Plotting of mne.combine_evoked output based on AverageTFR input (mutiple channels) only shows TFR output of the first channel

💬 Support & Discussions
1
  • MNE-Python version: 0.22.0
  • operating system: Windows 10

Hi everyone,

Please refer to the following figures:

2013_Combine1
2013_Combine11872×1152 258 KB

2013_Single
2013_Single1872×1152 214 KB

I’d like to draw attention to the middle plots (red-blue freq plot) in each of these two figures, which were created using the following codes extract:

ch_name = ['MEG2013','MEG2022','MEG2023'] #For top plot
ch_name = ['MEG2013'] #For bottom plot

epochsFIF_condtion1.pick_channels(ch_name)
epochsFIF_condtion2.pick_channels(ch_name)

tfr_epochs_1 = tfr_morlet(epochsFIF_condtion1, freqs,
                                          n_cycles=n_cycles, decim=decim, use_fft=False,
                                          return_itc=False, average=False)
                
tfr_epochs_2 = tfr_morlet(epochsFIF_condtion2, freqs,
                                          n_cycles=n_cycles, decim=decim, use_fft=False,
                                          return_itc=False, average=False)

this_tfr1 = tfr_epochs_1.average()
this_tfr2 = tfr_epochs_2.average()

## Appended as this is a loop function that creates tfr_morlet data out of multiple participants.
## all_power_sets are a collection of data of all participants
all_power_set.append(this_tfr1)
all_power_set2.append(this_tfr2)   

## Loop ends here after all participants' data have been extracted. Codes below this are only reached
## after the all_power_sets have been filled with all participants' data.
##---------------------------------

## Create averaged data across all participants
tfr_condition_1 = mne.grand_average(all_power_set, interpolate_bads=True, drop_bads=True)
tfr_condition_2 = mne.grand_average(all_power_set2, interpolate_bads=True, drop_bads=True)

tfr_contrast = mne.combine_evoked([tfr_condition_1, tfr_condition_2],
                                      weights=[1, -1])
tfr_contrast.plot(axes=ax2, vmin=vmin, vmax=vmax)

The top middle plot was created based on information from 3 sensors whereas the bottom middle plot was created based on information from 1 sensor. However, both plots look exactly the same. Changing the position of sensors in ch_name makes no difference as it always plot the sensor with the smallest number (ie. 2013). I was expecting the top plot to show the averaged time-freq activity across the three sensors but this is not happening.

Am I missing a step somewhere to get the plotting done right? Or is this some kind of bug?

2

you can compare numerical values with numpy and do the numerics looking at tfr_condition_1.data etc.
Basically the data attributes of the TFR objects.

but looking at the plots they look similar but not exactly the same.

Alex

3

Hi @agramfort,

Thanks for the reply! Still a little puzzled. I have now tried a different comparison using 1 vs. 20 gradiometers but the TFR output is looking the same. I believe that is unlikely that the TFR plot for 20 grads would be anywhere near the single grad.

compare
compare1241×1276 103 KB

Here above I have the TFR output from 1 grad and 20 grad, which looks too similar to me. The two tfr_contrast objects look correct in terms of the length of data corresponding to the number of sensors, as reflected in the console output. The data is right, but tfr_contrast.plot() for the 20 grads seem to be plotting the TFR of only the first array in the list of tfr_contrast object. Shouldn’t it be attempting to plot the average TFR of all 20 grads instead?

4

can you replicate using one dataset from mne? maybe the mne somato dataset?

so we can investigate.

Alex

5

There are several things that don’t work as expected; could the behavior you’re seeing be the issue described here?

(Or the related issues BUG: TFR.plot misc issues · Issue #8013 · mne-tools/mne-python · GitHub and https://github.com/mne-tools/mne-python/pull/9150)

6

Hey @cbrnr,

Yes, it seems like the same exact issue.

My tfr_contrast is an AverageTFR as well. So I guess the fix is still underway?

@agramfort I can try to replicate using the somato data but on a brief glance of the somato data file, it seems like there is only one event in the entire data file. Gonna go grab dinner for now. I can try to play around later in the evening by adding an additional baseline event to maybe make a simple working example closest to my combine_evoked setup, if required?

7

@AaronAng can you check out PR #9150 and see if this fixes the problem for you?

8

@cbrnr Tried replacing the tfr.py file with the following commit mne-python/tfr.py at tfr_misc_issues · eort/mne-python · GitHub but I’m getting an error. This was generated using the codes with the somato data example at the bottom of this reply.

error
error851×755 26 KB

@agramfort the issue can be replicated using the somato example. Note that plotting picks='all', is no different from plotting the first channel in picks=[0].

import os.path as op

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.time_frequency import tfr_morlet, psd_multitaper, psd_welch
from mne.datasets import somato

from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt')

data_path = somato.data_path()
subject = '01'
task = 'somato'
raw_fname = op.join(data_path, 'sub-{}'.format(subject), 'meg',
                    'sub-{}_task-{}_meg.fif'.format(subject, task))

# Setup for reading the raw data
raw = mne.io.read_raw_fif(raw_fname)
events = mne.find_events(raw, stim_channel='STI 014')

# picks MEG gradiometers
picks = mne.pick_types(raw.info, meg='grad', eeg=False, eog=True, stim=False)

# Construct Epochs
event_id, tmin, tmax = 1, -1., 3.
baseline = (None, 0)
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=baseline, reject=dict(grad=4000e-13, eog=350e-6),
                    preload=True)

epochs.resample(200., npad='auto')  # resample to reduce computation time

# define frequencies of interest (log-spaced)
freqs = np.logspace(*np.log10([6, 35]), num=8)
n_cycles = freqs / 2.  # different number of cycle per frequency
power, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, use_fft=True,
                        return_itc=True, decim=3, n_jobs=1)

## Able to plot specific sensor data
power.plot(picks=[82], baseline=(-0.5, 0), vmin=-0.6, vmax=0.6, mode='logratio', title=power.ch_names[82])
## But once you set picks='all', .plot() only plots the first channel in the list of channels.
power.plot(picks='all', baseline=(-0.5, 0), vmin=-0.6, vmax=0.6, mode='logratio', title='all')
## Notice that the first channel's activity is a replica of the 'all' plot.
power.plot(picks=[0], baseline=(-0.5, 0), vmin=-0.6, vmax=0.6, mode='logratio', title=power.ch_names[0])
## Even selecting a reduced list of sensors lead to the same error of plotting only the first sensor
power.plot(picks=['MEG 0113', 'MEG 0112', 'MEG 0122', 'MEG 0123', 'MEG 0132'], baseline=(-0.5, 0), vmin=-0.6, vmax=0.6, mode='logratio', title='sub-selected sensors')



mne.viz.tight_layout()
plt.show()
9

you’re answer is in

https://github.com/mne-tools/mne-python/blob/main/mne/time_frequency/tfr.py#L1383

@larsoner wrote this apparently.

can you open an issue on github maybe?

Alex

1 Like