Averaging n number of trials within Epoch object

đź’¬ Support & Discussions
,
1

Hi there,

I’m using macOS Catalina 10.15.7 and MNE 0.23.4

I am trying to average every 5 trials for each of the two conditions I have in my data after preprocessing to subsequently use in MVPA (based off Grootswagers et al., 2017) .

I have “AUD_Left” and “AUD_Right” conditions, each with N instances after equalising epochs.

The appearance AUD_Left and AUD_Right is random throughout the session.

I would like to average every 5 instances within a condition, so that I end up with N Ă· 5 instances of each condition after averaging.

Currently;

epochs_from_file

Number of events	190
Events	AUD_Left: 95
AUD_Right: 95
Time range	-0.500 – 0.998 sec
Baseline	-0.500 – -0.200 sec


X = epochs_from_file.get_data()

X.shape

(190, 64, 768)

y = epochs_from_file.events
y

array([216, 216, 216, 203, 203, 216, 203, 216, 216, 203, 216, 203, 203,
       203, 216, 203, 216, 203, 203, 203, 203, 203, 216, 216, 216, 216,
       216, 203, 216, 216, 216, 203, 203, 203, 216, 203, 203, 216, 203,
       216, 216, 216, 203, 203, 216, 203, 203, 216, 203, 216, 203, 203,
       203, 203, 216, 203, 216, 216, 216, 216, 203, 216, 216, 216, 216,
       203, 203, 203, 203, 216, 216, 203, 216, 203, 216, 216, 216, 216,
       203, 203, 216, 203, 216, 203, 203, 216, 216, 203, 203, 216, 216,
       203, 203, 203, 216, 203, 216, 203, 216, 203, 216, 216, 203, 216,
       216, 216, 203, 203, 203, 203, 216, 216, 203, 216, 216, 203, 203,
       203, 216, 216, 216, 216, 203, 203, 203, 216, 203, 216, 216, 216,
       203, 203, 203, 203, 203, 216, 203, 203, 216, 203, 216, 203, 216,
       216, 203, 216, 216, 216, 203, 203, 216, 203, 216, 203, 203, 203,
       203, 203, 216, 203, 216, 203, 203, 216, 203, 216, 216, 203, 203,
       216, 216, 216, 216, 216, 216, 216, 216, 216, 203, 203, 216, 203,
       216, 203, 203, 203, 203, 203, 216, 216], dtype=int32)

I’d like to average (not sliding window) every 5 instances of event 216 and 203 separately, so that afterwards;

X.shape

(38, 64, 768)

I hope that is clear and I’m sorry if this is very simple, I haven’t been able to find any solution so far!

Cheers,
Yohan

1 Like
2

Hello @YohanWards and welcome to the forum!

You can use NumPy slicing to achieve this: simply index the epochs via [start:stop:step]. Here’s a reproducible example:

# %%
import mne

sample_dir = mne.datasets.sample.data_path()
sample_fname = sample_dir / 'MEG' / 'sample' / 'sample_audvis_raw.fif'

raw = mne.io.read_raw_fif(sample_fname, preload=True)
raw.pick_types(eeg=True, stim=True)
raw.filter(l_freq=1, h_freq=40)

events = mne.find_events(raw, stim_channel='STI 014')
event_id = {
    'auditory/left': 1,
    'auditory/right': 2,
    'visual/left': 3,
    'visual/right': 4,
    'face': 5,
    'buttonpress': 32
}

epochs = mne.Epochs(
    raw, events=events, event_id=event_id,
    tmin=-0.2, tmax=0.5, baseline=(None, 0),
    preload=True
)

idx_start = 0  # first epoch
idx_stop = -1  # last epoch
idx_step = 5   # step width

epochs_auditory = epochs['auditory'][idx_start:idx_stop:idx_step]
epochs_visual = epochs['visual'][idx_start:idx_stop:idx_step]

evoked_auditory = epochs_auditory.average()
evoked_visual = epochs_visual.average()

mne.viz.plot_compare_evokeds({
    'every 5th auditory': evoked_auditory,
    'every 5th visual': evoked_visual,
})

image

Best wishes,
Richard

3

Hi Richard,

Thanks for that. Correct me if I’m wrong though, but that will take every 5th trial, and then average them all together to form 1 evoked object?

What I’m after is to average the first 5 epochs together, then the next 5, and so on, for each condition. So that in the end, I have a total of 19 epochs (that are each made up of the average of 5 epochs) for each condition that I can subsequently use in a MVPA.

The main issue that I’m having is doing this averaging within the epochs object, so that I can still call epochs_from_file to use on the MVPA as per your decoding tutorial (thank you so much for that tutorial by the way!). If you have any idea on how to implement this I’d greatly appreciate it.

Cheers
Yohan

4

you’ll have to do it manually with numpy and create an EpochsArray from the averaged trials.

Alex

5

You could try creating evoked objects for every 5 epochs, will this work for your use case?

# %%
import mne

sample_dir = mne.datasets.sample.data_path()
sample_fname = sample_dir / 'MEG' / 'sample' / 'sample_audvis_raw.fif'

raw = mne.io.read_raw_fif(sample_fname, preload=True)
raw.pick_types(eeg=True, stim=True)

events = mne.find_events(raw, stim_channel='STI 014')
epochs = mne.Epochs(raw, events=events, preload=True)

num_epochs_to_avg = 5
idx_from = list(
    range(0, len(epochs), num_epochs_to_avg)
)
idx_to = list(
    range(num_epochs_to_avg, len(epochs), num_epochs_to_avg)
)
idx_from_to = list(zip(idx_from, idx_to))

evokeds = []
for idx_from, idx_to in idx_from_to:
    evoked = epochs[idx_from:idx_to].average()
    evokeds.append(evoked)

evokeds[0].plot_joint()
evokeds[1].plot_joint()

image
image

2 Likes
6

Thank you Richard and Alexandre for the advice!

I think that code does want I want in terms of trial averaging.

However, the issue I’m facing now is using the evokeds object in my MVPA pipeline. Because it’s considered a list object I can’t extract an array like I have done below:


#Define variables to be used in cross validation 

X = epochs_from_file.get_data()
y = epochs_from_file.events[:,2]

#Define size of sliding window in samples (remember this is after downsampling) 
#So for example here "25" equals 25 samples at the new downsampled rate of 512 samples/s. Thus, 25 samples = 48.83ms)
winsize = 25

#Create "smoothedX" array to use to create sliding window

sX = np.empty(X.shape)

#Loop through channels/epochs and smooth the data by "winsize"

for ch_idx in range(X.shape[1]):
    for tr_idx in range(X.shape[0]):
        sX[tr_idx,ch_idx,:] = np.convolve(X[tr_idx,ch_idx,:], np.ones(winsize)/winsize, mode = 'same')
    
#Classifier pipeline.
clf = make_pipeline(
    # An MNE scaler that correctly handles different channel types –
    Scaler(epochs_from_file.info),
    # NumPy array reshaping thanks to the MNE vectorizer
    Vectorizer(),
    # And, finally, the actual classifier.
    SVC(kernel='linear'))

# Run cross-validation.
# Note that we're using MNE's cross_val_multiscore() here we simply pass the number of desired CV splits,
# and MNE will automatically do the rest for us.


scores = np.empty((5,sX.shape[2]))
for sample_idx in range(sX.shape[2]):
    scores[:,sample_idx] = (cross_val_multiscore(clf, sX[:,:,sample_idx], y, scoring='roc_auc'))

Is there any way I can turn this evokeds (which is the result of evoked.append) back into an epoch object, or an array of shape [38, 64, 768] so that I can call it in this pipeline?

7

You need to extract the data of each evoked as a NunPy array and then use all these arrays to create a new one. Something like this (untested, I’m on the phone):

data = np.array(
    [e.data for e in evokeds]
)