I calculated the power spectrum for multiple conditions, across multiple subjects, using the following code.

spectrum = epochs_cleaned.compute_psd(fmin=2, fmax=40., n_fft = 500, n_overlap = int(0.50 * 500), method = 'welch', n_jobs=-1)
# average across epochs first
mean_spectrum = spectrum.average()

I would like now to average the mean_spectrums calculated as below, across participants, and then plot each averaged condition as I was doing at the single subject / single condition level using

mean_spectrum.plot_topomap(cmap = 'RdBu_r')

How can I achieve this? As in, how can I aggregate different mean_spectrums and visualise the results?

Assuming that the conditionality of the data is the same across subjects (same number of channels, and same frequency axis) I see two ways around it:

(Probably discouraged) Get the “spectrum” or “mean_spectrum” (not sure if one of the two options would not work; I’ll assume that the 2nd works and use it in the example below) as you do for each subject and aggregate the results like

all_sub_data = []
all_sub_data.append([mean_spectrum.get_data() for sub in subjects])
all_sub_data = np.array(all_sub_data)
# average across subjects
all_sub_average = np.mean(all_sub_data, axis=0)

Then use one subject and replace its data in-place with the across-subjects averaged data mean_spectrum._data = all_sub_average.

The other way around it would be to get the data from each subject before computing the PSD, creating an EpochsArray object with the data from all subjects, then using the functions as in your example.