annotate_amplitude for a specific time window

Hello everyone,

for my EEG data I want to artifact reject based on an amplitude threshold (-100 to +100 ”V) and a joint probability threshold (SD = 5), but only for the first second after stimulus onset. My trials are 9 seconds long, and of course I can just epoch the data to the first second and run the annotation on it. But after this procedure I need my 9 sec epochs back. I am not sure how to tell python to integrate the corrected first second to the rest 8 seconds.

Thank you!

Franzi

You could:

  1. epoch your data in the first second
  2. run annotate_amplitude
  3. save the epoch indices that are considered bad
  4. epoch your data with your 9 seconds
  5. apply the indices you have saved before to reject the “bad” epochs

annotate_amplitude() only works with continuous (“raw”) data.

1 Like

right, thanks for the reminder Richard.

I think in that case I’d recommend using .get_data, numpy, and custom code.

Something like:

time_of_interest = epochs.time_as_index([0, 1])  # from 0 to 1 seconds
data = epochs.get_data()  # (n_epochs, n_channels, n_times)
for epoch in data:
    # epoch is a numpy array (n_channels, n_times)
    epoch_window = epoch[:, time_of_interest[0]:time_of_interest[1]]
    # epoch_window is the first second of your data
    ...  # and so on

Probably should use epochs.time_as_index() here to get the correct time indices!

Thank you for your help
I tried to run the following code

time_of_interest = data.time_as_index([0, 1])  # from 0 to 1 seconds
data1 = data.get_data()  # (n_epochs, n_channels, n_times)
for epoch in data1:
    # epoch is a numpy array (n_channels, n_times)
    epoch_window = data1[:, time_of_interest[0]:time_of_interest[1]]
    
reject_criteria = dict(eeg=100e-6)       
flat_criteria = dict(
                     eeg=-100e-6)

data_annot = mne.Epochs(epoch_window, events, tmin=0, tmax=1, reject_tmax=0,
                    reject=reject_criteria, flat=flat_criteria,
                    reject_by_annotation=False, preload=True)
data_annot.plot_drop_log()

It says undefined “events”. But when I do define the events it says "‘EpochsFIF’ object has no attribute ‘annotations’

time_of_interest = data.time_as_index([0, 1])  # from 0 to 1 seconds
data1 = data.get_data()  # (n_epochs, n_channels, n_times)
for epoch in data1:
    # epoch is a numpy array (n_channels, n_times)
    epoch_window = data1[:, time_of_interest[0]:time_of_interest[1]]
    
all_events, all_event_id = mne.events_from_annotations(data)    
event_id = slice(all_event_id, 'Stimulus/S ')
events, _ = mne.events_from_annotations(data, event_id=event_id)

reject_criteria = dict(eeg=100e-6)       
flat_criteria = dict(
                     eeg=-100e-6)

data_annot = mne.Epochs(epoch_window, events, tmin=0, tmax=1, reject_tmax=0,
                    reject=reject_criteria, flat=flat_criteria,
                    reject_by_annotation=False, preload=True)
data_annot.plot_drop_log()

Any idea how to handle this problem?

If I read your code correctly, you pass data to events_from_annotations(); however, this function expects a Raw and not an Epochs object.

Since you already managed to create epochs, there must be events somewhere??! You should just re-use them.

Best,

Richard

Yes I defined the events before running the ICA
I already tried to include the code here:

#define functions
def slice(sourcedict, string):
    newdict = {}
    for key in sourcedict.keys():
        if key.startswith(string):
            newdict[key] = sourcedict[key]
    return newdict

def flux_eegEpochs (raw_load):
    
    reject_criteria = dict(eeg=100e-6)       
  flat_criteria = dict(
                     eeg=-100e-6)

    all_events, all_event_id = mne.events_from_annotations(raw_load)
    print(all_event_id)
    print(all_events[:5])
   
    #Reload events chooseing only the ones that start with 'Stimulus/S '
    event_id = slice(all_event_id, 'Stimulus/S ')
    events, _ = mne.events_from_annotations(raw_load, event_id=event_id)

    # Epoch data
    epochs = mne.Epochs(raw_load, 
                    events, 
                    event_id=event_id, 
                    tmin=-1, 
                    tmax=8,
                    detrend=1,
                    reject=reject_criteria, 
                    flat=flat_criteria,
                    reject_by_annotation=False,
                    #baseline=(-0.2,0),
                    proj=True,
                    preload=True)

    len(epochs.events)
    
    return epochs

    epochs.plot_drop_log()

Doesnt work unfortunatly

Since your event names all start with Stimulus/, we can abuse MNE’s event name grouping capability: you can trivially select all Stimulus epochs via

epochs_stim = epochs['Stimulus']

No need to slice etc.

Does that help solve your problem?

Unfortunatly not. It does help to get rid of the “event” error but then it says "Flat value must be a number >= 0, not “-0.0001”, if I change it to “0” what I actually dont want because I want a -100 flat and +100 peak, it gives me the following error: “No appropriate channels found for the given picks (array(, dtype=int32))”

 epochs = mne.Epochs(raw_load, 
                    events, 
                    event_id=event_id, 
                    tmin=-1, 
                    tmax=8,
                    detrend=1,
                    reject=reject_criteria, 
                    flat=flat_criteria,
                    reject_by_annotation=False,
                    #baseline=(-0.2,0),
                    proj=True,
                    preload=True)

    len(epochs.events)
    
    return epochs

    epochs.plot_drop_log()

def identify_ica(raw):
    """ Use ICA to reject artifacts
    """
    # Perform ICA
    ica = mne.preprocessing.ICA(
            n_components=20, # Number of components to return
            max_pca_components=None, # Don't reduce dimensionality too much
            random_state=0,
            max_iter=800,
            verbose='INFO')
    ica.fit(raw, reject_by_annotation=True)
    
    # Plot ICA results
    
    ica.plot_components(inst=raw) # Scalp topographies - Click for more info
    ica.plot_sources(raw) # Time-courses - click on the ones to exclude
    
    return ica

raw_load=mne.io.read_raw_fif('C:\\Users\Admin\Desktop\Projekt_Franziska\Studie 3\Analyse\EEG\BIDS_mast\sub-56\ses-Expo\eeg/sub-56_ses-Expo_run01_data_interpolated.fif', verbose=False)
mne_data_folder ='C:\\Users\Admin\Desktop\Projekt_Franziska\Studie 3\Analyse\EEG\BIDS_mast\sub-56\ses-Expo\eeg/'

print(raw_load.annotations)
print(len(raw_load.annotations))
print(set(raw_load.annotations.duration))
print(set(raw_load.annotations.description))
print(raw_load.annotations.onset[0])

dataepoch=flux_eegEpochs (raw_load)
print(dataepoch.info)

filt_data = dataepoch.copy()
filt_data.load_data().filter(l_freq=1., h_freq=None)

raw_ica=identify_ica(filt_data);

Hello @FranziMNE, it’s become difficult for me to figure out which line exactly is causing the error you’re reporting. Could you please try to reduce your code snippet to the relevant lines only, and also share the entire traceback when you post error messages? This should make it easier for others to understand where things are going wrong.

Best wishes,
Richard

Got it, no epochs were left when I set the reject_criteria to the stringent threshold I used, so no epochs were left in the ICA, resulting in that error.
So I think everything can be traced back to the reject_criteria.

For the critera: Why cant I set the flat criteria to -100 microVolt? I actually want to have those criteria: amplitude threshold (-100 to +100 ”V) and a joint probability threshold (SD = 5)

The rejection thresholds are specified as so-called “peak-to-peak amplitude”, PTP.

This is basically the signal change from the lowest to the highest point. Therefore, the value can never be negative.

What you’re asking for is, to my knowledge, not supported. But I’m not sure if this is actually a problem — I believe the concept of PTP is actually better suited for artifact rejection, as it’s independent of signal baseline and DC offset.

Best wishes,
Richard

Alright! :slight_smile: The backround ist just, that for my current paper, reviewers asked to redo our EEG analysis since the preprocessing pipeline differs between ERP and TFR analysis in the paper. This is due the fact that my co-author ran the ERP analysis with matlab and I did the TFR analysis in python. I now set up the ERP analysis in python and try to build up the pipeline my co-author used. He used the mentioned artifaction correction criteria. This is why I wanted to redo them in python.
Thank you very much for your help again! :slight_smile:

Best, Franzi

I just don’t think there’s a trivial way to do it
 @agramfort may have an idea?

I think you mis-understand what flat means here. both reject and flat must be positive, and reflect absolute values of signal amplitudes. If you want to reject epochs that have values above 100 microvolts and also below -100 microvolts, then all you need to specify is reject=100 (actually, 100e-6, since MNE-Python stores data in Volts, not microvolts). flat is used to reject epochs where the signal stays below a given absolute value for the entire epoch (as an indication that something went wrong with the connection / the active electrode’s amplifier isn’t working / whatever).

Objection :wink:

It’s not signal amplitudes, it’s the peak-to-peak amplitudes! That is, the difference between the highest and the lowest value.

But your suggestion does do the trick for the scenario pointed out here, I believe, yes! So my previous statement:

is actually incorrect.

But I just wanted to clarify the nomenclature here a little.

ah, of course you are right. The main point stands (that flat is not appropriate here) but reject will not be sufficient by itself. I wonder if we should allow reject to take a callable?

Actually 
 no, this is not the case. Say if you specify reject={'eeg': 100e-6} and within a given epoch, the minimum signal amplitude is 40 ”V and the maximum is 120 ”V, the absolute amplitude exceeds 100 ”V but the PTP does not, so the epoch is retained.

We cannot “fix” this by specifying an additional flat parameter either, because this one also only looks at PTP.

I believe this would be the best approach that grants us greatest flexibility without forcing us to introduce a new parameter. @FranziMNE could then just use a function like:

def custom_reject(single_epoch_data):
    if (
        single_epoch_data().min() < -100e-6 or
        single_epoch_data().max() > +100e-6
    ):
        return False  # reject
    else:
        return True  # keep

Hi :slight_smile: Thanks for the idea

# Load data (after ICA)
subject = 'sub-52'
session = 'ses-Expo'
bids_root ='C:\\Users\\Admin\\Desktop\\Projekt_Franziska\Studie 3\\Analyse\\EEG\\BIDS_mast'
mne_data_folder=bids_root + '/' + subject + '/' +  session + '/' +  'eeg/' 
data=mne.read_epochs(mne_data_folder + subject  +  '_' + session + '_run01_' + 'data_ica.fif', verbose=False)


def custom_reject(data):
    if (
        data().min() < -100e-6 or
        data().max() > +100e-6
    ):
        return False  # reject
    else:
        return True  # keep

data.plot_drop_log()

It sadly doesnt reject anything, even when changing the criteria to very unreasonable low values