How to use mne.combine_evoked for subtraction?

[anapesq]

Hello all,

I can’t quite figure out how to use ‘mne.combine_evoked’ to subtract two conditions.
I want to subtract CON-INC. And when I test this, Iam getting the same exact result independently of the order of the weights [1, -1] vs [-1, 1].

• CON-INC

evoked_diff = mne.combine_evoked([CON, INC], weights=[1, -1])

• INC-CON

evoked_diff = mne.combine_evoked([CON, INC], weights=[-1, 1])

I am getting the same exact output from both options. Does anyone know where I am going wrong?

For background, here are the code lines and plots:

CONGRUENT

CON_average
Out[106]: <Evoked | 'Grand average (n = 25)' (average, N=25), 0 – 1 sec, baseline off, 204 ch, ~6.6 MB>

fig1 = CON_average.plot_topomap(ch_type='grad', time_unit='s',times=times_plt, colorbar=True,  vlim=(0,12),  cmap=cmap_v, nrows=2, ncols='auto',  contours=3)

INCONGRUENT

INC_average
Out[107]: <Evoked | 'Grand average (n = 25)' (average, N=25), 0 – 1 sec, baseline off, 204 ch, ~6.6 MB>

fig2 = INC_average.plot_topomap(ch_type='grad', time_unit='s',times=times_plt, colorbar=True,  vlim=(0,12) , cmap=cmap_v, nrows=2, ncols='auto',  contours=3)

CON-INC

evoked_diff_CONINC = mne.combine_evoked([CON_average, INC_average], weights=[1, -1])
Out[108]: <Evoked | 'Grand average (n = 25) - Grand average (n = 25)' (average, N=12.5), 0 – 1 sec, baseline off, 204 ch, ~6.6 MB>

fig3 = evoked_diff_CONINC.plot_topomap(ch_type='grad', time_unit='s',times=times_plt, colorbar=True,  vlim=(0,16),cmap=cmap_v, nrows=2, ncols='auto',  contours=3)

INC-CON

evoked_diff_INCCON = mne.combine_evoked([CON_average, INC_average], weights=[-1, 1])
Out[109]: <Evoked | '-Grand average (n = 25) + Grand average (n = 25)' (average, N=12.5), 0 – 1 sec, baseline off, 204 ch, ~6.6 MB>

fig1 = evoked_diff_INCCON.plot_topomap(ch_type='grad', time_unit='s',times=times_plt, colorbar=True,  vlim=(0,16),cmap=cmap_v, nrows=2, ncols='auto',  contours=3)

Many thanks!

Ana

[dasdiptyajit (Diptyajit Das)]

@anapesq Hi Ana, thanks for posting.
‘mne.combine_evoked’ API is used to subtract two evoked conditions. Here is the example code snippet, where aud_evoked and vis_evoked are two evoked conditions:
evoked_diff = mne.combine_evoked([aud_evoked, vis_evoked], weights=[1, -1])

If you reverse the the weight [-1, 1], you would get the same values/amplitudes as the output (i.e, as a part of the subtraction method) but with a sign reversal. Therefore, topography sign would reverse w.r.t the weights [1, -1]. You can easily test this with an eeg data set.

In your case, I would immediately sense that CON-INC vs INC-CON topographies should have a different sign polarity. But I can see from the figures that its not the case!

What kind of MEG system is it? Axial-gradiometers?? What is the output for: evoked_diff_CONINC.data Vs evoked_diff_INCCON.data . Do you see a sign reversal in the values?

best,

[drammock (Dan McCloy)]

On quick look, this seems like a bug in our plotting code to me. Do you have time to open a new bug report and include a reproducible example? Ideally you’d replicate it with one of the MNE sample datasets, if that can’t be done then we’ll need (download links for) the two files CON_average and INC_average. Please also include a link to this thread in the issue report.

[drammock (Dan McCloy)]

Meanwhile one thing you should try is using the defaults for vlim and cmap. I don’t think that’s the problem, but it doesn’t take long to check.

[anapesq]

Thanks for your reply @drammock.
I tried the vlim and cmap defaults and it didn’t solve the issue (picts below).
I have never submitted a new bug report, but will try to do that.

CON-INC

INC-CON

[drammock (Dan McCloy)]

@anapesq

I tried the vlim and cmap defaults and it didn’t solve the issue (picts below).

thanks for confirming.

[anapesq]

Thank you for your input @dasdiptyajit.

Good point! There is a sign reversal in the data. Would this point to a plotting bug as suggested by @drammock ?

The recording system is a MEGIN Neuromag (planar gradiometers)

evoked_diff_CONINC.data
Out[119]: 
array([[-2.98749794e-14, -3.56315137e-15,  2.74642740e-14, ...,
        -1.24641613e-13, -1.36985347e-13, -1.42341196e-13],
       [-9.85727203e-14, -1.96738810e-14, -1.62964724e-14, ...,

       

evoked_diff_INCCON.data
Out[120]: 
array([[ 2.98749794e-14,  3.56315137e-15, -2.74642740e-14, ...,
         1.24641613e-13,  1.36985347e-13,  1.42341196e-13],
       [ 9.85727203e-14,  1.96738810e-14,  1.62964724e-14, ...,

[dasdiptyajit (Diptyajit Das)]

If I try to recall my physics knowledge: Planner gradiometers compute the gradient of the magnetic field change (i.e., a pair of coils at each sensor location) that is tangential to the head surface. Therefore, I would mainly expect the positive values above the source topography in the sensor level. So, it doesn’t represent the actual source topography as it’s for magnetometers/axial gradiometers. So, if I imply the theory, I would expect CON-INC vs INC-CON topographies to be the same for planner gradiometers. @agramfort, will you correct me if I am wrong in this case! Thanks.

[richard (Richard Höchenberger)]

https://mne.tools/stable/generated/mne.Evoked.html#mne.Evoked.plot_topomap

ch_type ‘mag’ | ‘grad’ | ‘planar1’ | ‘planar2’ | ‘eeg’ | None

The channel type to plot. For 'grad', the gradiometers are collected in pairs and the RMS for each pair is plotted. If None the first available channel type from order shown above is used. Defaults to None.

The root-mean-squared value is always positive. It’s a deliberate choice for gradiometer visualization in MNE.

Best wishes,
Richard

[dasdiptyajit (Diptyajit Das)]

That makes more sense. Initially I thought those plots are from axial gradiometers but then its rather an issue with cmap parameter @anapesq- please use the defaults. Thanks for clarifying @richard.

[anapesq]

@dasdiptyajit Unfortunately, using the default cmap and vlim did not solve the issue. Please see plots for that in my previous reply.

[richard (Richard Höchenberger)]

As I said above, what you see is the expected behavior for gradiometer topomaps. If you have other channels like magnetometers, try plotting these and you should see a sign reversal.

Best wishes,
Richard

[anapesq]

Hi @richard, Thank you. Will do that!

[dasdiptyajit (Diptyajit Das)]

What you actually like to see (I assume a dipolar topomap that changes over time showing both positive and negative values) isn’t possible with gradiometers maps as richard mentioned. But if you have magnetometers along with gradiometers, you can simply remap/project the gradiometers data into magnetometers to reconstruct the dipolar topography.

Here is a nice example to follow:
https://mne.tools/dev/auto_examples/preprocessing/virtual_evoked.html

best
Dip