How do you correctly run multiple sensor data with permutation_cluster_test?

I have referred to the following guides on permutation_cluster_test:

https://mne.tools/stable/auto_tutorials/stats-sensor-space/plot_stats_cluster_time_frequency.html#sphx-glr-auto-tutorials-stats-sensor-space-plot-stats-cluster-time-frequency-py
https://mne.tools/stable/auto_tutorials/stats-source-space/plot_stats_cluster_time_frequency_repeated_measures_anova.html?highlight=permutation%20cluster%20test

They have been really helpful. However, the examples are based on the data off one single sensor and I would like to look at the average activation between a gradiometer pair, and possibly even a cluster of gradiometer pairs. As I am a novice to this statistical process, I need advice on the right way to go about doing this.

Below is an except of the codes from one of the guides that I believe is key:

ch_name = 'MEG 1332'

# Load condition 1
reject = dict(grad=4000e-13, eog=150e-6)
event_id = 1
epochs_condition_1 = mne.Epochs(raw, events, event_id, tmin, tmax,
                                picks=picks, baseline=(None, 0),
                                reject=reject, preload=True)
epochs_condition_1.pick_channels(ch_name)

# Load condition 2
event_id = 2
epochs_condition_2 = mne.Epochs(raw, events, event_id, tmin, tmax,
                                picks=picks, baseline=(None, 0),
                                reject=reject, preload=True)
epochs_condition_2.pick_channels(ch_name)    

decim = 2
freqs = np.arange(7, 30, 3)  # define frequencies of interest
n_cycles = 1.5

tfr_epochs_1 = tfr_morlet(epochs_condition_1, freqs,
                          n_cycles=n_cycles, decim=decim,
                          return_itc=False, average=False)

tfr_epochs_2 = tfr_morlet(epochs_condition_2, freqs,
                          n_cycles=n_cycles, decim=decim,
                          return_itc=False, average=False)

tfr_epochs_1.apply_baseline(mode='ratio', baseline=(None, 0))
tfr_epochs_2.apply_baseline(mode='ratio', baseline=(None, 0))

epochs_power_1 = tfr_epochs_1.data[:, 0, :, :]  # only 1 channel as 3D matrix
epochs_power_2 = tfr_epochs_2.data[:, 0, :, :]  # only 1 channel as 3D matrix

threshold = 6.0
T_obs, clusters, cluster_p_values, H0 = \
    permutation_cluster_test([epochs_power_1, epochs_power_2], out_type='mask',
                             n_permutations=100, threshold=threshold, tail=0)

If I would like epochs_power_1 and epochs_power_2 to each be a tfr output of a list of sensors, is it simply as easy as replacing ch_name = 'MEG 1332' to a list of sensors that I am interested in? Or is there more that I need to do? (ie. in epochs_power_1 = tfr_epochs_1.data[:, 0, :, :] # only 1 channel as 3D matrix, the comment made it seems like in order to read the data of multiple sensors, this might not be the right array output to do so?)

As far as I understand, epochs_condition_1.pick_channels(ch_name), when ch_name is a list of sensors, the proceeding tfr_epochs_1 output will be a tfr_morlet output that consists of the activity of the average activity from the list of sensors. Thus, epochs_power_1 = tfr_epochs_1.data[:, 0, :, :] would then be the single 3D matrix that contains the averaged activity. Is this right?

there is no clear “neighbor” structure between sensors that not from the same type. A gradiometer and a magnetometer at the same location can pick very different effects.

I fear I would advise to use one sensor type.

Alex

I see. What about selecting multiple sensors within the same sensor type then? (eg. a list [cluster] of gradiometers)? Assuming I’m interested in ‘MEG 1332’ and other gradiometers adjacent to 1332, do I just simply swap ch_name = 'MEG 1332' to ch_name = ['MEG 1332', 'MEG 1331', 'MEG 1342', 'MEG 1343', .....]?

I’ve tried doing that and the codes seem to be working fine - a final statistical plot was generated. I just wanted to be sure that the codes provided from the MNE guides are also capable of handling generation of statistics for averaged multiple sensors correctly.

1 Like