Hi! I am working with motor imagery dataset using mne, trying to calculate the ERDS ratio; I have used mne script from the website for ERDS for single participant, however, I am getting outliers (values are too high not sure why), I want to calculate the ratio for multiple participants, so trying to code the formulae by scratch i.e ERD/S = [(Pevent - Pref)/Pref ]* 100. Now, if I have to calculate Power from an epoch; for a specific channel, lets say Cz, how do I do that manually(I only know about the psd() function in mne)? Has anyone got a script for manual calculations!?
Or if I use the example listed in MNE, how do I remove the outliers and perform multiple subject average together? Thanks
I attach results of one person analysis, this looks a bit strange:
I am not exactly sure I understand what you are trying to achieve but here as some thoughts.
Regarding multiple subjects. I am not sure whether you need an across-subjects average baseline. If not, just loop the same code over subjects. If yes, I think you need to do it manually. Loop the code over subjects extracting the data with something like this_sub_data = epochs.get_data(), and in the final step average over all subjects.
Outlier rejection. I suppose you have some trials (and/or channels) with “outlier values”. This really depends on what you consider to be an outlier. A quick and dirty way would be for you to take a look into various epochs of the data, find and annotate bad segments of the data. An alternative, unsupervised approach would be to use the autoreject package which does that for you (but I would suggest you spend some time looking into what is rejected and why).
Manually computing power. Your formula is fine. Now regarding how to access the relevant data you would need something of this sort
# single subject baseline
baseline = epochs.copy().crop(-1.0, 0.0).pick_channels(['Cz'])
# Hilbert transform (specify bandwidth)
baseline = baseline.filter(<high_pass>, <low_pass>)
baseline = baseline.apply_hilbert()
# Hilbert transform power
Pref = baseline.get_data()
Pref = Pref ** 2
# same for Pevent with different time limits
event = epochs.copy().crop(0.0, 7.0).pick_channels(['Cz'])
...
@sotpapad Many thanks. This was really helpful. Actually I am trying to find an average ERD/ERS ratio w.r.to all the subjects.
I tried this, this is really helpful. I am stuck at the shape. I took Pref (-1,3.0) and Pevent (-1,9.). Now I have two different powers with different shapes.
Pref(40,1,2001) and size 80040 and , Pevent(40,1,5001) and size 200040.
How do I do the subtraction? Pref - Pevent? shall I reshape Pref as Pevent by adding series of 0s in Pref? or any other way? I am not confident with if this is the right way for substraction?
Since power evolves over time, I think that what most people would do is take the average during the baseline period and subtract it from every time point in the task period. In other words, you need Pref to be of dimensions (40,1). This way subtraction should work (and you can easily reshape it if you need to).
Of course, this may not be ideal especially if you do not consider your baseline to be somewhat stationary. In that case maybe you would prefer first computing a moving & weigthed average within x time points (e.g. every 100 ms) that would best capture local changes in power and then a median or something.
Hi @sotpapad Many thanks. This actually worked well for me. Event had a dimension of (1,5001) and Pref as (1,) . {I took overall average of trials}.
Many Thanks. Just to confirm, and be sure, please do correct me if I am wrong, when I took Hilbert transform, I had the data in real and imaginary part, when I apply Pref = Pref **2 , I did get real and imaginary that is the magnitude and phase (x+iy), now when I apply ERDS formula I expected % as I am looking at ratio. Before applying the ERDS formula do I just use the magnitude?
The envelope/magnitude would be abs(Pref) which is equivalent to using apply_hilbert(envelope='True'). Pref ** 2 corresponds to power spectral density (aka squared magnitude). And finally your formula gives you ERDS.
@sotpapad Many thanks . I tried it and getting the following result by plotting erd/ers manually: is this normal? I filter 8-32 Hz, remove bad epochs, and apply all that we discussed before and plot it. this motor imagery is average of few subjects at 0-9 sec window, event shown and happening 4 - 7. reference baseline is 0-3 sec.
Is this normal? Difficult to say. Maybe. I’ve data like this before, but that doesn’t really mean anything. It looks like it’s more noisy than what you would expect, but roughly counting some cycles I’d say it’s in the beta band alright. And, at least, there seems to be a ERD-like trend in your plot.
A couple of remarks that may help.
Remove epochs on the raw signals, then filter.
If you are using the autoreject package you will have different results depending on whether you provide an Epochs object that contains only channel Cz or multiple channels.
You probably need to discard the very first and last timepoints from your analysis as you obviously have some edge effects there.
I would start by plotting each subject separately, possibly with the original example code. Then proceed to the grand-average.
@sotpapad I tried removing bad epochs first and then filter and reduced time points from 2 to 8 second, I am taking baseline of 2.0 to 3.5 as Prest now and Pevent is 2.0 to 8.0 . I apply hillbert transform to get erd/ers and used envelope to get average. I get the following graph for channel C3:
I do not know what you are comparing those figures against, what data this is etc. However, the figure doesn’t look inherently wrong. It could be the case that this subject is (or the signals from those channels for this subject are) not that great in differentiating “left”.
There is really not much I can do to help you at this point. Sorry.
One final thing; I am unsure of the way you have implemented baseline correction. At some point you mentioned having a single baseline across trials and subjects. I would do it per trial and per subject.
