mass univariate analysis on individual electrodes and time points


  • Python version: 3.10.2
  • MNE version: 0.24.1
  • OS: Linux Manjaro 21.2.4

Hi everyone,

I looked for a similar question but couldn’t find it in the forum (apologies if I missed it). I am brand new to Python and MNE, so I’m possibly biting more than I can chew… any help would be greatly appreciated!

I have a preprocessed EEG dataset with 64 electrodes and two conditions, and I would like to identify ERP components in a data-driven way. Specifically, I want to average the stimulus-locked ERPs of both conditions and run permutation t-test on this grand-averaged dataset, in order to identify statistically significant activity at each electrode and time point. This would result in an unbiased localization of brain activity that should correspond to well-known ERP components (e.g., N1), based on the significant time window and scalp topography. In a second phase, I plan to use the electrode ROI and time window to test condition differences.

I can already do this in MATLAB with the Mass Univariate ERP Toolbox, i.e., calculate one-sample t-tests vs. 0 at each time point and electrode and control the familywise error rate using a t_max approach. According to this tutorial, the permutation_t_test function in mne should accomplish a similar goal.

Unfortunately, I am not able to retain the electrode dimension: the permutation_t_test function seems to automatically average activity across electrodes. Below a (hopefully) reproducible example:

###########################
# import libraries
###########################

import random
import numpy as np
import glob
import re

import mne
from mne.stats import permutation_t_test

###########################
# set parameters
###########################

project_seed = 999 # RNG seed
random.seed(project_seed) # set seed to ensure computational reproducibility

preproc_path = '/path/to/data' # directory with preprocessed files

filenames = glob.glob(preproc_path + '/*.fif') # list of .fif files in directory

conditions = ['cond1', 'cond2'] # condition names

###########################
# load data & compute evoked responses
# (this example only with dataset from one participant)
###########################

# load data, select scalp electrodes only
epochs = mne.read_epochs(filenames[0], preload = True).pick_types(eeg = True, exclude = ['M1', 'M2', 'VEOG', 'HEOG'])

# equalize epochs across conditions (drop epochs from condition with more data)
epochs_equal, dropped_epochs = epochs.equalize_event_counts()

###########################
# permutation t-test
###########################

# select scalp electrodes only
picks = mne.pick_types(epochs_equal.info, eeg = True, exclude = ['M1', 'M2', 'VEOG', 'HEOG'])

# separate condition-specific epochs
epochs_equal_cond1 = epochs_equal['cond1']
epochs_equal_cond2 = epochs_equal['cond2']

# extract all data and put them in one array
all_epochs = [epochs_equal_cond1.get_data(),
              epochs_equal_cond2.get_data()]

# define time window of interest
# in this example I don't need to run t-tests on all time points, because I'm only interested in the N1 (this also lowers the chance of Type II error)
t_min = 0
t_max = 0.3

# extract only time points of interest
temporal_mask = np.logical_and(t_min <= epochs_equal.times, epochs_equal.times <= t_max)

# average across conditions, then across epochs, then extract values only from time points of interest
all_epochs_perm = np.mean(np.mean(all_epochs, axis = 0), axis = 0)[:, temporal_mask]

# verify that the array dimensions are as expected (64 channels x 154 time points)
np.shape(all_epochs_perm)

# number of permutations for t-test
n_perm = 5000

# t-test (at last!)
T0, p_values, H0 = permutation_t_test(
    all_epochs_perm,
    n_permutations = n_perm,
    tail = 0,
    n_jobs = 1,
    seed = project_seed)

The expected output is (for both T0 and p_values) an array of dimensions 64 x 154 (channels x time points), one value per electrode x time point.

The actual output is (for both T0 and p_values) an array of dimensions 154 x 1 (time points), one value per time point :sob:

Am I using the wrong function? If not, is there a reason why permutation_t_test does not retain the electrode dimension? Is there a hack that would allow me to do that?

Thanks!

Hi everyone, please ignore the question above. Thanks to my collaborators, we figured out that I was incorrectly averaging across trials. Apologies for wasting your time!