ICA and line noise can be unstable?

[kingjr]

I have been trying to remove the line noise of some ECoG data using ICA (I am trying to avoid bandstop filter here).

However, the results I get are pretty unstable, and I am not quite sure why. That is, fitting the ICA several times of the same data can give radically different results (very good or very bad).

• I think I checked that this was not a copy/memory problem.
• Setting the ICA(random_state) stabilizes the results (but is arbitrary, since different random_state give very different results).

Here is a summary of what I do [complete notebook here]:

• sources = 1/F random source signal
• mix these source signals together with a line noise (raw = dot(cov, vstack(sources, linenoise))
• filtered_raw = raw.copy().filter(55, 65) to filter around line noise (because of this ref)
• ICA().fit(filtered_raw)`
• compute PSD of the data in component space dot(ica.components, original_raw)
• find the component carying the line noise: bad
• remove it from the data: ica.apply(original_raw, exclude=[bad])
• check the PSD of this cleaned data

This works ~ well for low sampling frequency with or without decim>1, but gets unstable for high sampling frequency when we use a decim>1.

This could be due to aliasing, so I tried to shuffle the data just prior to fitting the ICA

times = np.arange(raw.n_times)
np.random.shuffle(times)
raw._data = raw._data[:, times]

the results are a bit better, but still pretty unstable. In the notebook, you can see that it works 50% of the time

Any idea what I do wrong?

[agramfort]

I have been trying to remove the line noise of some ECoG data using ICA (I am trying to avoid bandstop filter here).

However, the results I get are pretty unstable, and I am not quite sure why. That is, fitting the ICA several times of the same data can give radically different results (very good or very bad).

welcome to the world of ICA :) it's a non-convex problem so sometimes you end up in a good local minimum sometimes you don't ...

besides if the spectral diversity of the source is what makes them different temporal ICA is probably not ideal ... but this would bring you far...

[mmagnuski]

Are you using method= ‘extended-infomax’? If I remember correctly this method should be better in catching line noise components.

[mmagnuski]

And just by the way: do you think implementing eeglab's cleanline in mne would be a worthwhile investment?

[jona-sassenhagen]

I like cleanline, although when I was using it, it had a really high crash propensity.  It has a bunch of advantages over Ica. 

[kingjr]

@mmagnuski extended-infomax seems more stable indeed, thanks!

Note that when using decim>1, we need to shuffle time.

Just had a quick look at cleanline, seems like a more robust approach in principle. Anyone willing to implement it in python?...

[jona-sassenhagen]

Not me :) Extended infomax was especially developed to deal with components like line noise which are non Gaussian.  At this point, my vote would be for deprecating standard infomax. 

[larsoner]

my vote would be for deprecating standard infomax

+1 from me, or at least actually using extended as the default as it sounds like we advertise it that way. Okay with you @dengemann?

Just had a quick look at cleanline, seems like a more robust approach in principle

Can you describe briefly what it does?

Have you tried "spectrum_fit" mode?

http://mne-tools.github.io/stable/generated/mne.io.RawFIF.html#mne.io.RawFIF.notch_filter

[jona-sassenhagen]

@Eric89GXL I'm too lazy to check right now, but in theory, it should always be clear if you're using extended or not because you call infomax via method='infomax' or method='extended-infomax'.

But it might be preferable to have both link to extended-infomax and only allow turning off extended manually.

[dengemann]

Okay with you @dengemann?

+0 ... mabye someone needs it for research questions. No pressure to deprecate it.

[jona-sassenhagen]

You can still set it via fit_params right?

[larsoner]

@dengemann even without deprecation of non-extended mode, do you agree extended should be on by default? It sounds like the default is non-extended even though extended is supposed to be the default based on the docs.

[jona-sassenhagen]

Just had a quick look at cleanline, seems like a more robust approach in principle

Can you describe briefly what it does?

It fits sinusoids at the expected line noise (and harmonics) to each channel and subtracts the best.

[dengemann]

do you agree extended should be on by default?

no, our default is fastica and it's the ICA method preferred by me and Alex.

[larsoner]

It fits sinusoids at the expected line noise (and harmonics) to each channel and subtracts the best.

The best what?

If it's doing an inner product and subtracting out the component, it sounds like it should be less robust than using "spectrum_fit" (which uses multitaper methods). Also, filter_chpi does this sort of thing (fits individual sisusoids) in a time-varying way, so you could use that if you want -- we'd just need to allow a mode that actually ignores cHPI signals, or generalize it to a different function name for this purpose.

[dengemann]

well maybe it is ok to kind of change the order of the documentation and name infomax-extended first and suggest to use that if fastica is not used.

[jona-sassenhagen]

@Eric89GXL I was wrong, but see:

Many analysts automatically perform a notch filter at 60 Hz to remove line noise. Such notch filters often use a notch width of 10 Hz or larger, resulting in significant signal distortion in frequencies between 50 and 70 Hz. This distortion should not be an issue for analyses that immediately apply a low-pass filter to the signal, say at 40 Hz, but may preclude certain high-frequency studies. Some studies have also shown that non-causal low-pass filtering can significantly change ERP onsets (Vanrullen, 2011; Rousselet, 2012). In addition, Barnett and Seth (2011) have shown that filtering, particularly low pass filtering, can have a damaging impact on Granger causality and other connectivity computations.

Mitra and Pesaran (1999) suggest a multi-taper decomposition approach for identifying and removing line noise components while minimizing background signal distortion. We leverage routines from the cleanline EEGLAB plugin developed by Mullen (2012), which extends functionality from the open source Chronux toolbox (Mitra and Bokil, 2007). This method traverses the data using a short sliding window (4 s with a 1-s slide by default). The method transforms the data in each window into the frequency domain using a set of Slepian or multi-tapers with a predetermined spectral resolution. Such tapers are ideal for isolating spectral energy within frequency bands, even for short time windows. By default, we use a taper bandwidth (TBW) of 2 Hz in each 4-s sliding window (W), which contains N sample points. The Slepian tapers are created by calling dpss(N, TBW_W/2, TBW_W–1) where dpss is the discrete prolate spherical sequence function from the MATLAB Signal Processing Toolbox.

The cleanline method fits a frequency-domain regression model to estimate the amplitude and phase of a deterministic sinusoid of a specified frequency, embedded in locally white noise. This is an idealized model of sinusoidal line noise of unknown phase and amplitude. A Thompson F-test assesses whether the complex amplitude is significantly non-zero. If the amplitude is significant (p < 0.01), the method reconstructs the time-domain sinusoid for the line noise frequencies. The method stitches together results from successive overlapping windows by using a sigmoidal weighted average specified by a smoothing parameter tau (100 by default). Finally, the method subtracts this fitted signal from the data. We repeat this process (a maximum of 10 iterations by default) until the sinusoid amplitude for the specified frequencies is not significantly different from the background.

In practice, the exact line frequency is unknown and variable. Following Mitra and Pesaran, we apply the regression model across a range of frequencies centered on each candidate line-noise frequency and select the frequency that maximizes the Thompson F-statistic. The range is +/− fScanBandwidth, which is 2 Hz by default. The user must provide a rough estimation of the line frequencies or use the defaults of 60 Hz with harmonics that are multiples of 60 up to the Nyquist frequency. The advantage of this approach over notch filtering is that it removes only deterministic line components, while preserving “background” spectral energy. The sliding window estimation further allows for non-stationarities in the phase and amplitude of the line component. The PREP pipeline function encapsulates this denoising functionality in the function cleanLineNoise.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4471356/

[dgwakeman]

Warning, this may be dgwakeman noise, but

Most line noise is not a pure sinusoid, so it ends up being a different signal from the cHPI in that respect

[kingjr]

@Eric89GXL I don't clearly understand spectrum_fit, an example would welp.

@jona-sassenhagen With regard to default changes, I think it is more critical to a detail the tutorials and explain when to use each ICA method, rather than changing the ICA defaults. ICA can be used in different contexts and for different purposes.

Most line noise is not a pure sinusoid

+1, and it's definitely not stationary.

[larsoner]

I don't clearly understand spectrum_fit, an example would welp. Mitra and Pesaran (1999) suggest a multi-taper decomposition approach

This part at least is what "spectrum_fit" is meant to implement. So you could look at the ref @jona-sassenhagen mentioned to see the relationships. I guess you could extend it to the "cleanline" approach if you want.

[jona-sassenhagen]

@jona-sassenhagen With regard to default changes, I think it is more critical to have a adapt or write clear tutorials on when to use each ICA method, rather than changing the defaults of ICA as ICA can be used in different contexts and for different purposes.

I've told Denis a few times I'm really impressed with MNE infomax, but IMO it's a mistake to make it so easy to use non-extended infomax. It's basically outdated, as far as I can tell. You ran into this issue because you assumed method='infomax' was a reasonable choice, which it isn't.

So basically, here is your tutorial: use 'fastica' or 'extended-infomax', depending on if you have the time to fix infomax and if you want your components ordered by explained variance. Don't use 'infomax'.

[larsoner]

Most line noise is not a pure sinusoid, so it ends up being a different signal from the cHPI in that respect

That is a good point -- it is not stationary, plus it has lots of harmonics. FYI cHPI is also not stationary, though -- head movements act as amplitude modulators on the cHPI signals, thereby spreading activation. That's why a time-varying fit is used in filter_chpi (which does the same thing as MaxFilter for cHPI artifact removal).

[kingjr]

@Eric89GXL If I understand correctly (else excuse for the discussion noise), spectrum_fit is used to notch out sinusoidal frequencies in each channel separately, isn't it? Whereas the ICA is multivariate: it aims at finding the common component that capture the line noise. In theory, there's no need to filter line noise/sinusoids once you regressed out the line component from the data.

[jona-sassenhagen]

You're losing rank though, and it assumes perfect stationarity and a few other things.

It's good, but not perfect - if I was looking into the frequencies around line noise, I'd probably be concerned.

[larsoner]

... Whereas the ICA is multivariate ...

I haven't looked into these much, but that sounds right to me.

[kingjr]

You're losing rank though

You're loosing rank, but in EEG it makes sense. One rank has to be devoted to the line/common reference, no?

it assumes perfect stationarity and a few other things.

I don't think it assumes frequency stationarity. The line noise could be disrupted, as long as the noise component is common across electrodes, the ICA will capture it. In fact, the notebook I shuffle time (and thus frequency stationariy) to make ICA more robust.

[jona-sassenhagen]

I mean spatially.

And with other methods you don't necessarily lose rank.

[mmagnuski]

@kingjr Do you need to shuffle time? Eeglab Infomax shuffles time by default (I'm not sure how often) but in general the order of samples shouldn't matter so much as components are assumed to be stationary. With respect to cleanline vs ica: I've had files where ICA removed line noise without much problem, but in some cases it was helpless. But I've had problems with cleanline too - I usually resort to cleanline when ICA doesn't help.

@Eric89GXL So it seems cleanline is very close to spectrum_fit. Cleanline is performed in time windows and the estimated parameters are then smoothed lineary or nonlineary (there's a parameter for that) across windows. This can be problematic at times btw.

[kingjr]

I mean spatially.

Ok. Is it common to have spatially varying line noise?

And with other methods you don't necessarily lose rank.

That makes me wonder whether we shouldn't actually keep bad ECoG (e.g. in white matter) to preserve a high rank (similar for loose EEG?).

Do you need to shuffle time?

If I use the decim of ICA, yes (see here), else you introduce aliasing before the ICA fit. But yes, ICA is blind to time.

[agramfort]

FYI the default EEGLAB solver is infomax

[jona-sassenhagen]

With extended 1 though. 

[larsoner]

Should we close this, or is there some code action to be taken? Maybe updating the ICA docs somewhere to mention that line noise can be challenging etc.?

[kingjr]

I'm closing it for now. Thank you all for your comments and help.