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?