Should average reference be recomputed after post-ICA channel interpolation?

Hi,

I am preprocessing BioSemi 128-channel EEG data in MNE-Python for ERP and sensor-level spectral analyses.

My pipeline is:

  1. Mark bad channels
  2. Average-reference excluding bads
  3. Fit/apply ICA on the good channels
  4. Interpolate bad channels using spherical splines

After interpolation, the data remain referenced to the average of the originally good channels, so is the mean across the reconstructed 128-channel montage is no longer exactly zero?

I therefore compared this pipeline with an additional average reference after interpolation. The overall spectral feautres remained highly correlated between pipelines, but baseline-corrected contrast magnitudes changed by approximately 8–30%, with larger changes in the participants with more/noisier bad channels.

Since the originally good channels already have a near-zero mean, the second average reference is mainly determined by the interpolated channels. Thus, its effect may depend on their number and spatial location.

My questions are:

  1. For ERP amplitudes and sensor-level spectral topographies, should the average reference be recomputed after interpolation (doing 2 CARs, one before ICA and one after interpolation), or should the reference to the originally good channels be retained?
  2. Does including interpolated channels in the final average undesirably over-weight their neighboring measured channels?
  3. Should the number and spatial distribution of interpolated channels be included as QC variables when using a post-interpolation average reference?

Thank you

Good points, and I think you might find the PREP pipeline helpful. After line noise removal, they basically detect and interpolate bad channels and then re-reference to the average of all good plus interpolated channels. After this basic preprocessing, you can continue with other processing steps like fitting ICA etc. If you want to give it a try, PyPREP is a Python package which implements exactly this pipeline (to be exact, there are some deliberate differences that should work even better than the original).

Hey Clemens, thanks for your response.

If I’m interested only in 0.5-45 Hz is actually needed to do line noise removal?, and also… since ICLabel wants the data average referenced, what would be the best approach; to do notch filter or to not?

Right now I’m doing two objects, one with (1, 100) filter to do ICA and then fitting the ICA into the object for analysis… as the following code shows:

raw.info['bads'] = bad_channels

raw.set_eeg_reference('average', ch_type='eeg') 

raw_ica = raw.copy().pick('eeg').filter(1, 100)

raw_analysis = raw.copy().pick('eeg').filter(0.5, 45)

n_components = 128 - len(bad_channels) - 1 # n_chan - bad_chan - 1 (avg_reference)

ica = ICA(n_components=n_components, max_iter='auto', random_state = 42, method='picard', fit_params=dict(extended=True, ortho=False))

ica.fit(raw_ica, picks='eeg')

ica.exclude = exclude

raw_clean = ica.apply(raw_analysis)

raw_clean.interpolate_bads(reset_bads=True)

But anyways, I’m still wondering if doing another average reference to the raw_clean final object is bad or not… I will give it a try to PyPrep.

You don’t have to remove line noise if you’re only interested in frequencies below 45 Hz. Regarding the second average reference, I would say it depends. Is important for your analysis that all participants have a theoretically comparable reference? Then I would do it. If not, no need. But PyPREP does give you a comparable reference, so if you want to go for comparable references, it is the even better approach compared to a second average reference in my opinion.

Thanks for your response I’m giving PyPrep a try now…

But also I wanted to ask you a question regarding ICA and ICLabel, would it be optimal to just do the ICA on the 0.5 - 45 copy object? Since when I’m passing the raw_ica object (1, 100) filtered with no line noise removal, my ICA decomposition gets a lot of components with line noise… so maybe I’m losing neural signal just because of not doing the ICA on the copy that I will really use… which approach would you recommend? maybe zapline in to the raw data and then separating the two objects for ICA?

Your response would be really helpful!

If you have strong line noise that spreads across many ICs, I would first remove it using a notch filter or a CleanLine/ZapLine-like approach (equivalent to mne.filter.notch_filter(method="spectrum_fit") – no need for another package unless MNE’s implementation doesn’t work well). Then fit the ICA on broad-band data (e.g., 1–100 Hz, or even without a low-pass filter using only a 1 Hz high-pass).

Using 0.5–45 Hz data for ICA is not ideal because muscle activity extends above 45 Hz, so restricting the frequency range can make separation more difficult. Also, ICLabel was trained on broad-band data, so components from narrow-band data may not match the distributions underlying its classes and the classifier will likely perform poorly.