ENH: combined gradiometers (plotting + stats)

[kingjr]

@dengemann @agramfort
Continuing on our email discussion... As I said, the pairs of gradiometers can be combined to get a 2D vector, and hence be treated as a single dependent variable.

In terms of plotting this could be translated into a 3D color plot where hues indicate the combined gradiometers orientations, and light to the norm of the vector.

For example (made from Matlab sorry...):

and in terms of topography

With such view
i) you're closer to what you record and (and actually appreciate the complexity of the collected data).
ii) you can tell when you have an inversion of polarity (although this is generally directly visible with magnetometers). For example, in http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0085791, we observed a sign reversal of the topographical pattern using MVPA. While the norm of combined gradiometers is obviously blind to such thing.

Regarding statistics, that would be nice to apply a single analysis using both the angle and the norm differ across two conditions, as
i) you would keep all information (unlike when you use the norm)
ii) you would reduce the number of tests (unlike a single statistical test per gradiometer). I haven't found such analysis yet. You could off course test this separately using circular analysis for the angle.

Philipp Berens suggested me to use bootstrapping for testing both at once, but on my computer that rapidly became unpractical using so many sensors and time points.

An interesting follow up question could then be related to the time frequency transform applied to such vector.

[kingjr]

Oh, and also, in MNE, the combined values are RMS, whereas in fieldtrip it's just the norm of the vector (http://fieldtrip.fcdonders.nl/reference/ft_combineplanar). I'm not sure which one is the most sensible...

[agramfort]

sounds like a good motivation with relevant usecases

you just have to port this code to Python, come up with a good API and open
the PR :)

looking forward to it !

are you around in June to come to the MNE sprint?

https://github.com/mne-tools/mne-python/wiki/Coding-sprints

On Tue, May 13, 2014 at 11:12 AM, kingjr [email protected] wrote:

@dengemann https://github.com/dengemann @agramforthttps://github.com/agramfort
Continuing on our email discussion... As I said, the pairs of gradiometers
can be combined to get a 2D vector, and hence be treated as a single
dependent variable.

In terms of plotting this could be translated into a 3D color plot where
hues indicate the combined gradiometers orientations, and light to the norm
of the vector.

For example (made from Matlab sorry...):

[image: example_time]https://cloud.githubusercontent.com/assets/4881164/2955175/fedda99c-da7c-11e3-8394-e8cd0adb3e97.png

and in terms of topography

[image: example_topo]https://cloud.githubusercontent.com/assets/4881164/2955189/2db8e178-da7d-11e3-9d1d-60114dbe8ea0.png

With such view
i) you're closer to what you record and (and actually appreciate the
complexity of the collected data).
ii) you can tell when you have an inversion of polarity (although this is
generally directly visible with magnetometers). For example, in
http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0085791,
we observed a sign reversal of the topographical pattern using MVPA. While
the norm of combined gradiometers is obviously blind to such thing.

Regarding statistics, that would be nice to apply a single analysis using
both the angle and the norm differ across two conditions, as
i) you would keep all information (unlike when you use the norm)
ii) you would reduce the number of tests (unlike a single statistical test
per gradiometer). I haven't found such analysis yet. You could off course
test this separately using circular analysis for the angle.

Philipp Berens suggested me to use bootstrapping for testing both at once,
but on my computer that rapidly became unpractical using so many sensors
and time points.

An interesting follow up question could then be related to the time
frequency transform applied to such vector.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/1280
.

[kingjr]

Yes. I am still learning the MNE/python structure (mitten coding style)... But I am decided to stop Matlab as fast as possible.
No I won't be here - somewhere in Canada... Next time..

[agramfort]

looking forward to see what you can contribute to the community :)

[dengemann]

@kingjr @agramfort I think I need this functionality sooner than later. Would be glad to put API issues first and later add different methods of combining grads.
I think the first thing we should add is a function that returns combined grads as a numpy array + a list of labels that indicate the origin of the channels. E.g. we would loose 2 channels for one bad gradiometer. We can later think about wrapping this in a meaningful info / fiff compatible structure.

WDYT?

Apart from that I'm wondering whether we should support merging grads at all stages. Also I'm not quitesure how to deal with analyses where you compute e.g. power spectral densities or TFRs. Would merging be recommended here, if so, at which stage?

to me the API of a designated Mixin would be as follows:

CombineChannelsMixin(object):
   """Mixin to combine channels"""
   def get_grads_merged(self, method='norm'):
       # returns numpy array 
       # XXX now
       pass

   def combine_channels(self, method='norm'):
       # returns instance with channels combined + info updated
       # XXX later
       pass 

[agramfort]

to me combined grads should be done at the evoked stage. We can had it
later to epochs if requested.

evoked_combined = evoked.combine_grads(method='norm', copy=True)

seems like a nice API. Then make everything private. You can take/refactor
the code from viz that does this.

the pain will be to update the viz code to work with evoked_combined.plot()

you'll need to create a new channel type. The new channel name can be
"MEGXX1-MEGXX2" ie the concat of the combined channels.

my 2c

[kingjr]

If I understand correctly, you want to output the norm of the combination of each pair of gradiometers for each trial/epoch? If that's the case, you'll have to be careful, because you'll square the within-trial noise too, and thus loose SNR. Furthermore, using the norm won't let you apply anymore time frequency transform on your signal.

This is why I suggested to remain in a complex space by default (grad1+1i*grad2). In such case, the norm will have to be retrieved after having first averaged data across trials: i.e. norm(mean(combined_data,axis="trials")). This step also allows you to later do time frequency analyses etc.

(If you go for that second option, I suggest that the new channels rather be called MEGXX1+MEGXX2.)

[agramfort]

agreed. norm on epochs is dangerous.

I vote for :

evoked_combined = evoked.combine_grads(method='norm', copy=True)

with method='norm' | 'angle' | 'complex'

what I fear is that having complex values in evoked.data will break a lot of code that expects evoked (IO, viz, etc...)

do we really need complex? can we just support norm + angle?

[dengemann]

@kingjr @agramfort let's focus on averages first. I like the proposed API but i would maybe start with some method that just returns the merged data instead of an MNE object. It would allow us to start somewhere and tackle all the issue mentioned as we go (e.g. naming, viz module, channel location information, fiff writing).
I would start with this method that returns a numpy array.

WDYT?

[agramfort]

Ok

[dengemann]

We could also start with a function. One method less to deprecate later ...

[adykstra]

I rather like @kingjr's suggestion of supporting complex values. @agramfort, couldn't there be a check in the modules expecting evoked that if type(evoked.data) is complex, do norm?

[dengemann]

@agramfort we can have complex values if we first go for a function.

[agramfort]

@agramfort https://github.com/agramfort we can have complex values if
we first go for a function.

ok

Andy we can go for complex but be prepared to change a bit of code
everywhere to avoid weird error messages

[agramfort]

@kingjr do you have time to make progress on this ? or clarify the way you want things to be done?

[kingjr]

I lost track of this issue. I'll try to rediscuss this with @dengemann next week if he's interested

[kingjr]

@agramfort This is apparently not on our priority list for now. The idea is simple, but the implementation can be complex, and would greatly benefit from virtual sensors which I don't master. So I'd rather postpone it.

[agramfort]

@kingjr feel free to reopen if you feel like it

[drammock]

There was a request for this feature on the gitter channel today; reopening

[finetard]

is there an implementation of this?

[drammock]

is there an implementation of this?

not yet. the closest we get is described here: https://mne.tools/stable/auto_examples/preprocessing/plot_virtual_evoked.html

[apesquita]

Hello all, Do you have new recommendations for combining planar gradiometers?
I would like to achieve something similar to this fieldtrip procedure: https://www.fieldtriptoolbox.org/reference/ft_combineplanar/
My ultimate goal is to visualize and analyze TFRs.
Thank you in advance for any pointers.

[jyeatman]

@drammock @larsoner - I wanted to circle back to see what the resolution was on this issue and what the current recommendation is for combining grads

[agramfort]

no progress here :(

while we do this internally for plotting we never give you access to some Raw, Evoked or Epochs with combined grads as it's technically a new ch_type (or we hack it as still a planar grad...) So it has consequences on plotting functions etc...

[enniin]

Hi, I'm also trying to analyse and visualise TFRs using grads. Just being able to get combined data as a numpy array sounded like a good idea to me! That is anyhow what e.g. mne.stats.permutation_cluster_test() requires.

[mscheltienne]

Adding a +1 as it was again requested recently on our site.