plot_compare_evokeds from mne.EvokedArray: confidence intervals and ylim

đź’¬ Support & Discussions
,
1

Platform: Windows-10-10.0.19041-SP0
Python: 3.7.6 (default, Jan 8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
mne: 0.23.0

Dear MNE users,

I am new with MNE and I am trying to plot ERPs data processed with EEGLAB (i.e., a matlab file with 128 rows -channels- and 513 columns -time points-). I have created a mne.EvokedArray object and could plotted it fine (see the attached figure) with the following code:

import mne
import scipy.io

samplesfile = scipy.io.loadmat("path/mat_file.mat")  #File read into dictionary
samples = samplesfile['data_mean2']  #Extract the numpy array in the dictionary

biosemi_montage = mne.channels.make_standard_montage('biosemi128')
length_channels = len(biosemi_montage.ch_names)

sampling_freq = 512
ch_types = ['eeg'] * 128
info = mne.create_info(ch_names=biosemi_montage.ch_names, ch_types=ch_types, sfreq=sampling_freq)

raw = mne.EvokedArray(samples, info, tmin=-.2, nave=160, kind='average', comment='Tone')
raw.set_montage(biosemi_montage)

mne.viz.plot_compare_evokeds(raw, picks=('C23'), colors=dict(Tone='red'), linestyles=dict(Tone='solid'))

I have only two questions I couldn´t solve:

  1. How can I plot CI of the ERP? Given that my data file contains data already averaged across subjects, I considered to use a callable object for CI. But when I tryed to add some random data with a callable object (see the example below), none shadow CI appeared around the ERP.
import numpy as np
def ci95(): return np.random.randint(500000000, size=(128, 513))

mne.viz.plot_compare_evokeds(raw, picks=('C23'), colors=dict(Tone='red'), linestyles=dict(Tone='solid'),
                            ci=callable(ci95))
  1. I don´t know why my data is scaled to those values in the Y axis (please, see the same attached figure). When plotting the topomaps, I could apply scalings=dict(eeg=1) and convert the Y lims of the colorbar to [-3, 3]. But I cannot figure it out how to scale it on the plot_compare_evokeds figure. I saw in the documentation that the parameter ylim makes

Y-axis limits for plots (after scaling has been applied).

But when I set this parameter in the plot_compare_evokeds as, for instance, ylim=dict(eeg=[-20, 20]), I got an error that says:

Image size of 510x54826907 pixels is too large. It must be less than 2^16 in each direction.

I am sure I am missing something, but I have read many forums and documentation and I couldn´t solve it.

Thanks in advance for any help on these issues.
Best,

Fer

figure_sample

2

First, a mathematical note: if all you have is the already-averaged data (that is all you show in your code sample, anyway), there is no way to know the confidence interval. Assuming, then, that your samplesfile structure has some other info in it like data_ci, then you would construct a function that conforms to what the docstring asks for:

If a callable, it must take a single array (n_observations Ă— n_times) as input and return upper and lower confidence margins (2 Ă— n_times).

In other words, your function should be something like:

def my_ci(array):
    ci = samplesfile['data_ci']  # or whatever your source of the CI info is
    assert ci.ndim == 2
    assert ci.shape[0] == 2
    assert ci.shape[1] == array.shape[-1]
    return ci

all those assert statements are not strictly necessary, but they enforce what is required of the my_ci function in order for it to work with plot_compare_evokeds. Finally, you would pass my_ci like this: plot_compare_evokeds(..., ci=my_ci) (in other words, don’t include () after the function name, and don’t wrap it with callable() either.

In MNE-Python, EEG data are stored in Volts but are usually converted to microVolts when plotting. If the values seem way off to you, perhaps your variable samples needs to be converted to Volts before you use it in mne.EvokedArray.

2 Likes
3

Dear Dan,

Thank you very much for your reply on my post. Regarding data scaling in the plot (i.e., question 2 above), I have solved it following your suggestion.

However, I am still a bit confused into how plotting CI in plot_compared_evokeds from a mne.EvokedArray object.

As you said, my object does not includes variability data, but I can easily generate this data in the same or another object. So, following your suggestion, I am not sure what should I use as input in:

Does it has to be just a similar array from that with averaged data across subjects (i.e., with the same size of N channels x N time points) with standard deviation or standard error data? Or does it has to be an array including only the specific channel I am trying to plot? Furthermore, I am not sure if it has to be just an array or a mne.EvokedArray( ..., kind='standard_error') object.

I have tried with some of this options above-mentioned but no shadowed CI appears in the plot. Thank you very much in advance for any help!

Best,

Fer

4

the ci argument works best in 2 specific use cases:

  1. you pass in a list of Evoked objects (say, one from each subject) and get an across-subjects confidence band
  2. you pass in Epochs.iter_evoked() (which converts each epoch into a separate Evoked object, sort of like each one is a one-item average), and you get an across-trials confidence band

If what you’re starting with is a single EvokedArray object and you don’t have access to the separate trials before they were averaged together then you have no way of generating a real confidence interval. If you do have the original data prior to averaging, then pass that to the plot_compare_evokeds function as in (1) or (2) above. If you don’t have the unaveraged data, but you do have pre-computed lower and upper confidence values for each timepoint, then you can shoe-horn that data into the plot in the way described in my previous answer (by creating a mock function that basically ignores the input array and returns the pre-computed (2 x n_times) array.

No. As mentioned in the docstring and in my previous answer, the shape of the CI should be (2 x n_times). In understanding why, it will help to remember that plot_compare_evokeds will always combine information across channels (see the combine parameter) and the CI is not really meant to show across-channel confidence bands. If you’re trying to plot a single channel, then your CI would be the average signal value in that channel ± the standard deviation (computed across replicates, which is presumably either across subjects or across trials). Or if you don’t want to plot standard deviation, use standard error or bootstrapped 95% CI or… take your pick.

1 Like
5

Dear Dan,

thank you very much again for your reply and all the explanation. I have now a clearer picture on how CI works in plot_evoked_compared. I am aware that my input object is just an array with data already averaged across subjects and trials, but I was trying to do exactly what you commented above:

However, and I still do not know why, I cannot make it work. I loaded an array of 2 x n_times (i.e., the same time points as the mean data) with upper and lower CI (also tried inverting lower and upper CIs values in the rows) as:

sample_CI_file = scipy.io.loadmat("path/data_CI_OneChannel.mat") 
sample_CI = sample_CI_file['data_CI_OneChannel']

And then called it with the function you suggested as:

def my_ci():
    ci = sample_CI['data_CI_OneChannel']
    assert ci.ndim == 2
    assert ci.shape[0] == 2
    assert ci.shape[1] == array.shape[-1]
    return ci

I tried many times in different ways but I couldn’t make it work (i.e., the code run fine but still no shadowed CIs are plotted).

Anyway, I will keep looking what is going on. Thank you again for all your help!

Best,
Fer

6

It turns out that the function has a bug in it, and it’s not drawing the CI if there’s only one condition being plotted. I’ll work on a fix for it.

2 Likes
7

Thank you very much, Dan!

Indeed I am trying to plot two conditions. I wrote here the minimum information possible for the problem I had and I was checking with one condition if it worked.

So, for two conditions, should I pass a list to ci as:

CI_A_file = scipy.io.loadmat("path/data_CI_OneChannel_condA.mat")
CI_A = CI_A_file['data_CI_OneChannel_condA.mat']

CI_B_file = scipy.io.loadmat("path/data_CI_OneChannel_condB.mat")
CI_B = CI_B_file['data_CI_OneChannel_condB.mat']

CI_conditions = [CI_A, CI_B]

def my_ci():
    ci = CI_conditions
    assert ci.ndim == 2
    assert ci.shape[0] == 2
    assert ci.shape[1] == array.shape[-1]
    return ci

I don´t see that is working either, but probably I am doing something wrong…
Please, apologize me for so many questions on this issue. Thank you once again.

Best,

Fer

8

The fix is in the works: Fix plot compare evokeds by drammock · Pull Request #9663 · mne-tools/mne-python · GitHub

This will still not work because you still are writing a function that doesn’t take any arguments. The first line MUST look like this:

def my_ci(array):  # or call the input "x" or "foo" or whatever you want

That is what it means when the documentation says

If a callable, it must take a single array (n_observations Ă— n_times) as input and return upper and lower confidence margins (2 Ă— n_times).

However, that won’t be enough to solve your problem (which it turns out is different than what you originally stated). Since you want to plot multiple conditions, that makes it rather tricky to hack together a function that will return the correct array of precomputed CI bounds for each condition that you’re trying to plot. Put another way, the rather unorthodoxed approach I showed in an earlier answer only really works for plotting single conditions.

Do you really not have access to the un-averaged data? It would be much easier to just pass the un-averaged data to plot_compare_evokeds and tell it what percentage of confidence interval you want (ci=0.95 or ci=0.68 or whatever)

1 Like
9

Dear Dan,

thank you once again for your reply. I am sorry for having posted an example (i.e., plotting CI for one condition) different of the real problem (i.e., to plot CI for multiple conditions). I thought (from reading the documentation) that plot_compared_evokeds could use a callable object for plotting CI either for a single condition or multiple conditions.

Anyway, I will try to do it from the unaveraged data or someway different. Thank you once again for your help.

Fer

10

no problem. You are right, the documentation did say that, and it was wrong (but the function has been fixed now, it will work to use a callable when you have a single condition because the pull request I linked to above has been merged). So if you’re able to use the current development version of MNE-Python (see here) then you can try it out. But since what you want to do is plot 2 conditions anyway, it may not be necessary. I still think your best bet is to pass in the unaggregated data, if you have access to it.