Creating adjacency matrix issue

Dear all,
Hope everyone is doing well!

I am trying to run a permutation clustering test, thanks to the last office hour, I managed to create a 3D data structure for it, but now I got an error message related to adjacency. Below you can find part of my code and the error message.

Thanks in advance for your help!

Cheers,
Gözem

data_path = '//'

file_names = fnmatch.filter(os.listdir(data_path), "*.fif")

for file_counter in file_names[0:len(file_names)]:
    data_fname = data_path + file_counter
    epoch = mne.read_epochs(data_fname, preload=True)

adjacency, ch_names = find_ch_adjacency(epoch.info, ch_type='eeg')

event_id = {'Stimulus/S 700': 700, # incorrect violation
                       'Stimulus/S 701': 701, # correct violation
                       'Stimulus/S 710': 710, # incorrect non-violation
                       'Stimulus/S 711': 711} # correct non-violation 

epoch.drop_channels(['HEL', 'HER', 'VED', 'ML', 'MR'])

epoch.equalize_event_counts(event_id)

adjacency, ch_names = find_ch_adjacency(epoch.info, ch_type='eeg')

evokeds = epoch.average(by_event_type=True)

n_subjects = 35
all_evokeds = [evokeds] * n_subjects

Xs = []
for k in range(len(evokeds)):
    X = np.array([e[k].data for e in all_evokeds])
    Xs.append(X)

The error message when I try to run spatio_temporal_cluster_test function:

adjacency (len 59) must be of the correct size, i.e. be equal to or evenly divide the number of tests (94459).
If adjacency was computed for a source space, try using the fwd[“src”] or inv[“src”] as some original source space vertices can be excluded during forward computation

  • MNE 1.3 Windows 11

There are some issues with your code sample, see below, with my comments in ALL CAPS:

data_path = '//'

file_names = fnmatch.filter(os.listdir(data_path), "*.fif")

for file_counter in file_names:  # NO NEED TO SAY [0:len(file_names)]
    data_fname = data_path + file_counter
    epoch = mne.read_epochs(data_fname, preload=True)
    #  THIS IS THE END OF YOUR FOR-LOOP. YOUR VARIABLE `epoch` WILL
    # ONLY CONTAIN THE DATA OF THE VERY LAST SUBJECT, BECAUSE THAT
    # VARIABLE IS GETTING OVERWRITTEN ON EACH ITERATION OF THE LOOP.
    # PROBABLY YOU SHOULD INDENT EVERYTHING THAT FOLLOWS.

# REMOVE THE NEXT LINE, IT'S REPEATED AGAIN LATER
adjacency, ch_names = find_ch_adjacency(epoch.info, ch_type='eeg')

event_id = {'Stimulus/S 700': 700, # incorrect violation
                       'Stimulus/S 701': 701, # correct violation
                       'Stimulus/S 710': 710, # incorrect non-violation
                       'Stimulus/S 711': 711} # correct non-violation 

epoch.drop_channels(['HEL', 'HER', 'VED', 'ML', 'MR'])

epoch.equalize_event_counts(event_id)

# MOVE THIS LINE OUTSIDE THE FOR-LOOP. FOR THE CLUSTERING TO WORK,
# YOU NEED THE SAME SET OF CHANNELS FOR EVERY SUBJECT, SO THE 
# ADJACENCY ONLY NEEDS TO BE CALCULATED FOR ONE OF THEM.
adjacency, ch_names = find_ch_adjacency(epoch.info, ch_type='eeg')

# THIS LINE SHOULD HAPPEN INSIDE THE FOR-LOOP, BUT THERE IS A PIECE
# MISSING: YOU NEED A CONTAINER (PROBABLY A `list`) CREATED BEFORE THE
# FOR-LOOP, AND YOU NEED TO ADD EACH SUBJECT'S `evokeds` TO THAT CONTAINER
# (IF IT'S A LIST, USE `all_evokeds.append(evokeds)` OR SIMILAR).
evokeds = epoch.average(by_event_type=True)

# THIS IS WHERE THE INDENTATION OF THE FOR-LOOP SHOULD STOP, PROBABLY



# THIS IS AN ODD CHOICE, YOU'RE COPYING ONE SUBJECT'S DATA AND TREATING IT
# AS THOUGH IT WERE 35 DIFFERENT SUBJECTS.
n_subjects = 35
all_evokeds = [evokeds] * n_subjects

Xs = []
# HERE YOU'LL WANT TO USE `range(len(all_evokeds))`
# HOWEVER, I THINK YOU'RE COLLAPSING ALL YOUR CONDITIONS TOGETHER,
# SO YOU'LL NEED TO TAKE A LOOK AT THAT TOO.
for k in range(len(evokeds)):
    X = np.array([e[k].data for e in all_evokeds])
    Xs.append(X)

Sorry that this doesn’t directly answer the question you asked… once you’ve fixed these things ping me and I’ll take another look.

1 Like

Hi Dan,

Many thanks, Dan for checking the code in detail! It helped a lot!
So, based on your feedback, I have rewritten it again. Do you think that makes sense?

event_id = {'Stimulus/S 700': 700, # incorrect violation
                   'Stimulus/S 701': 701, # correct violation
                   'Stimulus/S 710': 710, # incorrect non-violation
                   'Stimulus/S 711': 711} # correct non-violation 

all_evokeds_700 = list()
all_evokeds_701 = list()
all_evokeds_710 = list()
all_evokeds_711 = list()

for f in file_names:
    epoch = mne.read_epochs(data_path + f)
    epoch.equalize_event_counts(event_id)
    epoch.drop_channels(['HEL', 'HER', 'VED', 'MR'])
    adjacency, ch_names = find_ch_adjacency(epoch.info, ch_type='eeg')
    for e in event_id:
        if e == 'Stimulus/S 700':
            all_evokeds_700.append(epoch[e].average())
        if e == 'Stimulus/S 701':
            all_evokeds_701.append(epoch[e].average())
        if e == 'Stimulus/S 710':
            all_evokeds_710.append(epoch[e].average())    
        elif e == 'Stimulus/S 711':
            all_evokeds_711.append(epoch[e].average())   

grand_average_700 = mne.grand_average(all_evokeds_700)
grand_average_701 = mne.grand_average(all_evokeds_701)
grand_average_710 = mne.grand_average(all_evokeds_710)
grand_average_711 = mne.grand_average(all_evokeds_711)

ga_700_transpose = np.transpose(grand_average_700.data, (1,0))
ga_701_transpose = np.transpose(grand_average_701.data, (1,0))
ga_710_transpose = np.transpose(grand_average_710.data, (1,0))
ga_711_transpose = np.transpose(grand_average_711.data, (1,0))

X = test = [ga_700_transpose, ga_701_transpose, ga_710_transpose, ga_711_transpose]

ok, the code looks much better now. So let’s talk about the permutation test. Normally this is not done on grand averages, but rather you keep the separate subjects in there. Based on what you’ve shown so far, I’m guessing that this is a “within subjects” comparison, i.e., you’re interested in comparing “correct vs incorrect” or “correct violation vs correct non-violation”; all participants experienced all of these conditions.

If that is correct, then you want to do a one sample test (also called a paired differences test). This means you’ll construct the X array as a subtraction: something like [[subj1 correct violation], [subj2 correct violation], ... [subj_n correct violation]] *MINUS* [[subj1 correct nonviolation], [subj2 correct nonviolation], ... [subj_n correct nonviolation]]. In other words, you don’t want to call grand_average() before constructing X. Since I think you’re not doing any frequency analysis (right?) you can use mne.stats.spatio_temporal_cluster_1samp_test — MNE 1.3.1 documentation which means your X should have shape (n_subjects, n_times, n_channels). This means that you need to transpose the individual subject’s data instead of the grand-average data. Something like this:

# for the "correct violation versus correct nonviolation" comparison
X = (np.array([s.get_data().transpose() for s in all_evokeds_701]) - 
     np.array([s.get_data().transpose() for s in all_evokeds_711]))

The adjacency you created should work (it should have the right shape), and you can control clustering along the time axis by using the max_step parameter.

1 Like

That worked perfecty! Thanks a lot!