- MNE version: 1.0.3
- operating system: Windows 11
I am trying to plot the visualisation of the differences between the band frequency activity between two conditions. Referring to the attachments at the end of this post however, the display of the raw (normalized) activities seemed to be incongruent with one another. For instance, one of the few examples would be that the alpha plots are wildly different from one another, with the permutation plot showing negative t-values in areas where the raw activity (normalized) plot showed positive value. Logically speaking, if the raw activity is showing in one direction, it’s corresponding t-value should be in the same direction as well.
Does anyone have any idea why this might be the case?
Below is a snippet of the relevant codes directly pertaining to the question:
## Define Experiment Conditions
conditions = ['311','312']
Contrast1 = conditions[0]
Contrast2 = conditions[1]
epochs = mne.Epochs(
filt_raw, events, [int(item) for item in conditions], tmin, tmax, proj=False,
picks='eeg', baseline=baseline, preload=False, event_repeated='merge',
reject=reject, reject_by_annotation=True)
psd1, freqs = mne.time_frequency.psd_welch(epochs[Contrast1], fmin=1, fmax=40, reject_by_annotation=True)
psd2, _ = mne.time_frequency.psd_welch(epochs[Contrast2], fmin=1, fmax=40, reject_by_annotation=True)
bands = [(Delta[0], Delta[1], 'Delta'), (Theta[0], Theta[1], 'Theta'), (Alpha[0], Alpha[1], 'Alpha'), (Beta[0], Beta[1], 'Beta'), (Gamma[0], Gamma[1], 'Gamma')]
def plot_psd_topomap_diff_norm(psd1, psd2, bands, info, cmap):
psd_diff = psd1.mean(axis=0) - psd2.mean(axis=0)
fig, axes = plt.subplots(1, len(bands), figsize=(15, 5))
for i, (fmin, fmax, band_name) in enumerate(bands):
band_indices = np.where((freqs >= fmin) & (freqs <= fmax))[0]
band_power_diff = psd_diff[:, band_indices].mean(axis=1)
# Normalize the band power difference
band_power_diff_normalized = (band_power_diff - np.mean(band_power_diff)) / np.std(band_power_diff)
mne.viz.plot_topomap(band_power_diff_normalized, info, axes=axes[i], show=False, contours=5, cmap=cmap)
axes[i].set_title(f'Norm P(blue)-F(red) {band_name} ({fmin}-{fmax} Hz)')
plt.tight_layout()
plt.show()
# Plot the difference between conditions
plot_psd_topomap_diff_norm(psd1, psd2, bands, filt_raw.info, cmap)
def compute_band_power(psd, freqs, bands):
n_epochs, n_channels, n_freqs = psd.shape
band_powers = np.zeros((n_epochs, n_channels, len(bands)))
for i, (fmin, fmax, band_name) in enumerate(bands):
band_indices = np.where((freqs >= fmin) & (freqs <= fmax))[0]
if len(band_indices) == 0:
raise ValueError(f"No frequency indices found for band {band_name} ({fmin}-{fmax} Hz)")
band_power = psd[:, :, band_indices].mean(axis=2) # Mean over the frequency axis
band_powers[:, :, i] = band_power
return band_powers
# Compute band power for each condition
band_power1 = compute_band_power(psd1, freqs, bands)
band_power2 = compute_band_power(psd2, freqs, bands)
# Perform permutation test for each sensor and each band
p_values = np.zeros((band_power1.shape[1], len(bands)))
T_obs_all = np.zeros((band_power1.shape[1], len(bands)))
for ch in range(band_power1.shape[1]):
for band_idx, band in enumerate(bands):
data_ch_band = [band_power1[:, ch, band_idx], band_power2[:, ch, band_idx]]
# Perform the permutation cluster test
T_obs, clusters, cluster_p_values, H0 = permutation_cluster_test(
data_ch_band, n_permutations=1000, threshold=None, tail=0)
# Store the test statistics and p-values
T_obs_all[ch, band_idx] = T_obs[0] if len(T_obs) > 0 else 0
p_values[ch, band_idx] = cluster_p_values[0] if len(cluster_p_values) > 0 else 1
# Get sensor names
sensor_names = epochs.info['ch_names']
print('Sorted as: Delta Theta Alpha Beta Gamma')
print('-----------------------------------------------------')
# Print each row of p_values with the corresponding sensor name
for ch in range(len(sensor_names)):
formatted_p_values = " ".join(f"{p_val:.3f}" for p_val in p_values[ch])
print(f"Sensor: {sensor_names[ch]:<5} p-values: {formatted_p_values}")
def plot_significant_sensors(T_obs_all, p_values, bands, sensor_names, epochs, cmap):
for band_idx, (fmin, fmax, band_name) in enumerate(bands):
significant_sensors = np.where(p_values[:, band_idx] < 0.05)[0]
if len(significant_sensors) > 0:
data_to_plot = T_obs_all[:, band_idx]
# Create a mask for significant sensors
mask = np.zeros(data_to_plot.shape, dtype=bool)
mask[significant_sensors] = True
# Color code: red for condition 1 > condition 2, blue for condition 2 > condition 1
color_data = np.where(data_to_plot > 0, data_to_plot, -data_to_plot)
plt.figure()
sig_sensor_names = [sensor_names[s] for s in significant_sensors]
title = f'Significant sensors in {band_name} (t-values: R=Fail>Pass, B=Pass>Fail): {", ".join(sig_sensor_names)}'
plt.title(title)
# Plot the topomap
mne.viz.plot_topomap(color_data, epochs.info, mask=mask, cmap=cmap, show=True)
else:
print(f"No significant sensors found in {band_name} band.")
plot_significant_sensors(T_obs_all, p_values, bands, sensor_names, epochs, cmap)
