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?