MRG, BUG, DOC: derive GFP from average reference for EEG

[hoechenberger]

The reason for #8772 is simply that we always calculated the GFP incorrectly, i.e. without re-referencing sensor signals to average first. When doing that, I now correctly get a flat GFP (0 at all time points) for single-channel recordings:

This also means that since the introduction of this feature in 0.15 or so, our GFP plots haven't been correct in terms of absolute amplitude (which isn't too bad since it's difficult to judge the absolute amplitudes from the plot anyway)

Closes #8772, #8774

[hoechenberger]

Actually due to the nonlinearity in the term, even the relative changes in GFP were incorrect all along, if I'm not mistaken.

[hoechenberger]

I discussed this with @agramfort and therefore want to provide literature evidence why GFP must use the average-referenced signals.

TL;DR: GFP is based on average-referenced signal.

When using GFP, commonly two papers are cited:

• D. Lehmann, W. Skrandies, Reference-free identification of components of checkerboard-evoked multichannel potential fields, Electroencephalography and Clinical Neurophysiology, Volume 48, Issue 6, 1980, Pages 609-621, https://doi.org/10.1016/0013-4694(80)90419-8
• Dietrich Lehmann, Wolfgang Skrandies, Spatial analysis of evoked potentials in man—a review, Progress in Neurobiology,
  Volume 23, Issue 3, 1984, Pages 227-250, https://doi.org/10.1016/0301-0082(84)90003-0

In the 1980 paper, the term "global field power" does not yet occur; instead, it is called the "root-mean square", and they introduce the term "power":

In the 1984 paper, the term "global field power" shows up, albeit only for the "reference-free" calculation; they mention the perks of using the average-referenced signal in the paragraph below the formula:

In: Skrandies, W. Global field power and topographic similarity. Brain Topogr 3, 137–141 (1990), https://doi.org/10.1007/BF01128870, the GFP is explicitly named as being based on the average-referenced signal:

More recently, in a 2008 review tutorial by Murray MM, Brunet D, Michel CM. Topographic ERP analyses: a step-by-step tutorial review. Brain Topogr. 2008 Jun;20(4):249-64. doi: 10.1007/s10548-008-0054-5, the GFP is presented as average-reference based too:

Lastly, in a 2011 paper by Brunet et. al, GFP is also defined as average-reference-based (Brunet, Denis, Micah M. Murray, and Christoph M. Michel. "Spatiotemporal analysis of multichannel EEG: CARTOOL." Computational intelligence and neuroscience 2011, https://doi.org/10.1155/2011/813870)

[cbrnr]

What's the difference to #8774 (which I just merged)?

[hoechenberger]

#8774 was not supposed to be merged anymore ;) Because I figured it doesn't address the underlying problem: the incorrect calculation of the GFP, which is addressed in this PR

[hoechenberger]

Ok I reverted the merge of #8774

@agramfort Your call on this one :)

[hoechenberger]

Here a comparison between the approach currently used in master and the one implemented in this PR for EEG data:

import os
import matplotlib.pyplot as plt
import numpy as np

import mne


sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False)
events = mne.find_events(raw, stim_channel='STI 014')
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
              'visual/right': 4}
epochs = mne.Epochs(raw, events, tmin=-0.3, tmax=0.7, event_id=event_dict,
                    preload=False)
evoked = epochs['auditory/left'].average()
evoked.pick_types(meg=False, eeg=True)

fig, ax = plt.subplots(figsize=(12, 7))

ax.plot(np.sqrt((evoked.data ** 2).mean(axis=0)),
        label='master', ls='--')
ax.plot(evoked.data.std(axis=0, ddof=0),
        label='avg-ref', ls='--')

ax.set_xlabel('Sample', fontsize='large', fontweight='bold')
ax.set_ylabel('Global Field Power in V', fontsize='large', fontweight='bold')
ax.legend()

[agramfort]

to add some info to the debate. In https://neuroimage.usc.edu/brainstorm/Tutorials/Averaging I read: The green line represents the global field power (*GFP*), ie the sum of the square of all the sensors at each time point. This measure is sometimes used to identify transient or stable states in ERP/ERF. You can hide it with the display options menu *Extra > Show GFP*. can you check in brainstorm, fieldtrip or eeglab how it's done? Remark: this is no such things as common average reference on MEG as you don't need a reference :)

…

[hoechenberger]

The green line represents the global field power (GFP), ie the sum of the square of all the sensors at each time point.

But this would be something entirely different again! Unless there's a bug in their docs. They don't even take the root of the value, so it's actual power. GFP as outlined in all the publications I cited above is actually not a measure of power, as one can see by the resulting unit (µV, not µV^2).

can you check in rainstorm, fieldtrip or eeglab how it's done?

I don't remember a built-in function in EEGLAB; I always had to use a plugin / additional toolbox. On the mailing list, they recommend to simply calculate the standard deviation across channels, as I'm doing in this PR. This automatically uses the average-referenced signal.

Re Fieldtrip, Brainstorm I can do the research.

We could also discuss simply offering different "GFP" calculation methods.

[hoechenberger]

Remark: this is no such things as common average reference on MEG as you don't need a reference :)

Yep I'm aware of that, however you could still subtract the mean 😅

[hoechenberger]

@agramfort The FT docs are helpful regarding my previous-to-last comment (GFP not actually being a power):

FT_GLOBALMEANFIELD calculates global mean field amplitude or power of input data

Use as
[gmf] = ft_globalmeanfield(cfg, data)

The data should be organised in a structure as obtained from the
FT_TIMELOCKANALYSIS function. The configuration should be according to
FT_PREPROCESSING function. The configuration should be according to

cfg.method = string, determines whether the amplitude or power should be calculated (see below, default is 'amplitude', can be 'power')
cfg.channel = Nx1 cell-array with selection of channels (default = 'all'),
see FT_CHANNELSELECTION for details

This function calculates the global mean field power, or amplitude,
as described in:
Lehmann D, Skrandies W. Reference-free identification of components of
checkerboard-evoked multichannel potential fields. Electroencephalogr Clin
Neurophysiol. 1980 Jun;48(6):609-21. PubMed PMID: 6155251.

Please note that to calculate what is clasically referred to as Global
Mean Field Power, cfg.method must be 'amplitude'. The naming implies a
squared measure but this is not the case.

That aside, this is how it's calculated:

switch cfg.method
  case 'amplitude'
    dataout.avg = std(dataout.avg,1);
  case 'power'
    dataout.avg = std(dataout.avg,1).^2;
end

https://github.com/fieldtrip/fieldtrip/blob/3ba8125c1e809a8a9e6ed4f6f286b8fb3357d28a/ft_globalmeanfield.m#L95-L100

I'm not really sure what dataout.avg is exactly, but it's clear that for the amplitude GFP (the one I'm implementing in this PR, which is what is classically known as GFP), they simply calculate the standard deviation, i.e. use an average-referenced signal.

[hoechenberger]

@agramfort I'd propose to introduce an Evoked.gfp() method that gets parameters to control its behavior (use average-referenced signal or not, return power or amplitude).

[hoechenberger]

ERPLAB uses the standard deviation, too:

        MGFP_data   = std(bindata(:,:,j));

https://github.com/lucklab/erplab/blob/13ce4c384b52b472713c586e2454902fac57e818/functions/mgfperp.m#L77

[hoechenberger]

Updated tutorial:
https://24697-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/evoked/plot_eeg_erp.html#global-field-power-gfp

cc @drammock

[drammock]

some suggestions for making the tutorial a little tighter. Mostly nitpicks except the point about fixing ylim for the two plots that we're meant to compare

[hoechenberger]

I couldn't get viz/tests/test_evoked.py::test_plot_evoked_reconstruct to pass locally, and I'm not sure why. Would appreciate any help.

[larsoner]

@hoechenberger the values changed for EEG (rightly so) and so we can just adjust the limits in the test, will push a commit. Good to know that the test is sufficiently sensitive to changes in values :)

[hoechenberger]

will push a commit

Thank you, @larsoner!

[hoechenberger]

Tutorials look good:
https://24771-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/evoked/plot_20_visualize_evoked.html#plotting-signal-traces

https://24771-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/evoked/plot_eeg_erp.html#global-field-power-gfp

API docs look good:
https://24771-1301584-gh.circle-artifacts.com/0/dev/generated/mne.Evoked.html?highlight=evoked%20plot#mne.Evoked.plot

This is good to merge once green.

[hoechenberger]

Tutorials look good:
https://24771-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/evoked/plot_20_visualize_evoked.html#plotting-signal-traces

Whoops. No, It doesn't. The amplitude of the RMS is super big… the heck? scratches head

[hoechenberger]

Seems np.linalg.norm() uses the Frobenius norm:

This doesn't calculate the mean of the squared values before taking their square-root – right? Is there another norm we could use? Otherwise I'll just go back to the original approach, np.sqrt((D**2).mean(axis=0))

[larsoner]

you can use norm and divide by sqrt(n) or use the old formula, those two should be equivalent

[hoechenberger]

@larsoner I reverted to the old formula now because it's more intuitive to grasp when looking at the code. If you think using the approach based on the Frobenius norm would have strong advantages, I can switch to that too!

[hoechenberger]

All good now.
https://24776-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/evoked/plot_20_visualize_evoked.html#id1

[larsoner]

Thanks @hoechenberger !

[hoechenberger]

@larsoner Can you try to backport if you have a minute? Otherwise I can do it later today.

[larsoner]

done

[hoechenberger]

You're an angel <3 Or maybe a saint? who knows! Thank you!