Determine the significance of connectivity estimates against baseline connectivity

This example demonstrates how surrogate data can be generated to assess whether connectivity estimates are significantly greater than baseline.

# Author: Thomas S. Binns <t.s.binns@outlook.com>
# License: BSD (3-clause)
# sphinx_gallery_thumbnail_number = 3

import matplotlib.pyplot as plt
import mne
import numpy as np
from mne.datasets import somato

from mne_connectivity import make_surrogate_data, spectral_connectivity_epochs

Background

When performing connectivity analyses, we often want to know whether the results we observe reflect genuine interactions between signals. We can assess this by performing statistical tests between our connectivity estimates and a ‘baseline’ level of connectivity. However, due to factors such as background noise and sample size-dependent biases (see e.g. [1]), it is often not appropriate to treat 0 as this baseline. Therefore, we need a way to estimate the baseline level of connectivity.

One approach is to manipulate the original data in such a way that the covariance structure is destroyed, creating surrogate data. Connectivity estimates from the original and surrogate data can then be compared to determine whether the original data contains significant interactions.

Such surrogate data can be easily generated in MNE using the make_surrogate_data() function, which shuffles epoched data independently across channels [2] (see the Notes section of the function for more information). In this example, we will demonstrate how surrogate data can be created, and how you can use this to assess the statistical significance of your connectivity estimates.

Loading the data

We start by loading from the Somatosensory dataset, MEG data showing event-related activity in response to somatosensory stimuli. We construct epochs around these events in the time window [-1.5, 1.0] seconds.

# Load data
data_path = somato.data_path()
raw_fname = data_path / "sub-01" / "meg" / "sub-01_task-somato_meg.fif"
raw = mne.io.read_raw_fif(raw_fname)
events = mne.find_events(raw, stim_channel="STI 014")

# Pre-processing
raw.pick("grad").load_data()  # focus on gradiometers
raw.filter(1, 35)
raw, events = raw.resample(sfreq=100, events=events)  # reduce compute time

# Construct epochs around events
epochs = mne.Epochs(
    raw, events, event_id=1, tmin=-1.5, tmax=1.0, baseline=(-0.5, 0), preload=True
)
epochs = epochs[:30]  # select a subset of epochs to speed up computation

  0%|                                               | 0.00/611M [00:00<?, ?B/s]
  1%|▎                                     | 5.08M/611M [00:00<00:11, 50.8MB/s]
  2%|▋                                     | 10.2M/611M [00:00<00:12, 47.1MB/s]
  2%|▉                                     | 15.1M/611M [00:00<00:12, 48.2MB/s]
  3%|█▎                                    | 20.7M/611M [00:00<00:11, 51.3MB/s]
  4%|█▌                                    | 26.0M/611M [00:00<00:11, 51.7MB/s]
  6%|██                                    | 33.6M/611M [00:00<00:09, 59.9MB/s]
  7%|██▌                                   | 42.1M/611M [00:00<00:08, 68.1MB/s]
  8%|███▏                                  | 50.6M/611M [00:00<00:07, 73.4MB/s]
 10%|███▋                                  | 59.3M/611M [00:00<00:07, 77.8MB/s]
 11%|████▏                                 | 67.5M/611M [00:01<00:06, 79.1MB/s]
 12%|████▋                                 | 75.8M/611M [00:01<00:06, 80.0MB/s]
 14%|█████▎                                | 84.4M/611M [00:01<00:06, 81.9MB/s]
 15%|█████▊                                | 93.3M/611M [00:01<00:06, 84.2MB/s]
 17%|██████▌                                | 102M/611M [00:01<00:05, 85.2MB/s]
 18%|███████                                | 111M/611M [00:01<00:05, 85.7MB/s]
 20%|███████▋                               | 120M/611M [00:01<00:05, 87.0MB/s]
 21%|████████▏                              | 129M/611M [00:01<00:05, 87.3MB/s]
 22%|████████▊                              | 137M/611M [00:01<00:05, 85.6MB/s]
 24%|█████████▎                             | 146M/611M [00:01<00:05, 84.6MB/s]
 25%|█████████▉                             | 155M/611M [00:02<00:05, 85.5MB/s]
 27%|██████████▍                            | 163M/611M [00:02<00:05, 84.0MB/s]
 28%|██████████▉                            | 172M/611M [00:02<00:05, 85.3MB/s]
 30%|███████████▌                           | 181M/611M [00:02<00:05, 85.5MB/s]
 31%|████████████                           | 189M/611M [00:02<00:05, 82.7MB/s]
 32%|████████████▋                          | 198M/611M [00:02<00:04, 83.5MB/s]
 34%|█████████████▏                         | 206M/611M [00:02<00:04, 84.0MB/s]
 35%|█████████████▋                         | 215M/611M [00:02<00:04, 81.3MB/s]
 37%|██████████████▏                        | 223M/611M [00:02<00:04, 81.2MB/s]
 38%|██████████████▊                        | 231M/611M [00:02<00:04, 82.2MB/s]
 39%|███████████████▎                       | 240M/611M [00:03<00:04, 82.4MB/s]
 41%|███████████████▊                       | 248M/611M [00:03<00:04, 83.9MB/s]
 42%|████████████████▍                      | 257M/611M [00:03<00:04, 84.8MB/s]
 44%|████████████████▉                      | 266M/611M [00:03<00:04, 85.6MB/s]
 45%|█████████████████▌                     | 275M/611M [00:03<00:03, 86.2MB/s]
 46%|██████████████████                     | 284M/611M [00:03<00:03, 87.3MB/s]
 48%|██████████████████▋                    | 293M/611M [00:03<00:03, 88.2MB/s]
 49%|███████████████████▎                   | 302M/611M [00:03<00:03, 88.9MB/s]
 51%|███████████████████▊                   | 311M/611M [00:03<00:03, 89.2MB/s]
 52%|████████████████████▍                  | 320M/611M [00:03<00:03, 89.6MB/s]
 54%|████████████████████▉                  | 329M/611M [00:04<00:03, 89.8MB/s]
 55%|█████████████████████▌                 | 338M/611M [00:04<00:03, 89.9MB/s]
 57%|██████████████████████▏                | 347M/611M [00:04<00:02, 90.1MB/s]
 58%|██████████████████████▋                | 356M/611M [00:04<00:02, 90.3MB/s]
 60%|███████████████████████▎               | 365M/611M [00:04<00:02, 90.3MB/s]
 61%|███████████████████████▉               | 374M/611M [00:04<00:02, 90.1MB/s]
 63%|████████████████████████▍              | 383M/611M [00:04<00:02, 90.4MB/s]
 64%|█████████████████████████              | 392M/611M [00:04<00:02, 90.4MB/s]
 66%|█████████████████████████▋             | 401M/611M [00:04<00:02, 89.2MB/s]
 67%|██████████████████████████▏            | 410M/611M [00:04<00:02, 87.4MB/s]
 69%|██████████████████████████▊            | 419M/611M [00:05<00:02, 86.1MB/s]
 70%|███████████████████████████▎           | 427M/611M [00:05<00:02, 85.1MB/s]
 71%|███████████████████████████▊           | 436M/611M [00:05<00:02, 73.9MB/s]
 73%|████████████████████████████▎          | 444M/611M [00:05<00:02, 75.7MB/s]
 74%|████████████████████████████▉          | 453M/611M [00:05<00:02, 78.4MB/s]
 76%|█████████████████████████████▍         | 461M/611M [00:05<00:01, 81.0MB/s]
 77%|█████████████████████████████▉         | 470M/611M [00:05<00:02, 68.3MB/s]
 78%|██████████████████████████████▍        | 477M/611M [00:05<00:01, 68.1MB/s]
 79%|██████████████████████████████▉        | 484M/611M [00:06<00:01, 67.1MB/s]
 80%|███████████████████████████████▎       | 491M/611M [00:06<00:01, 67.8MB/s]
 82%|███████████████████████████████▊       | 498M/611M [00:06<00:01, 65.3MB/s]
 83%|████████████████████████████████▎      | 505M/611M [00:06<00:01, 66.7MB/s]
 84%|████████████████████████████████▊      | 514M/611M [00:06<00:01, 72.6MB/s]
 86%|█████████████████████████████████▎     | 522M/611M [00:06<00:01, 76.9MB/s]
 87%|█████████████████████████████████▉     | 531M/611M [00:06<00:00, 79.9MB/s]
 88%|██████████████████████████████████▍    | 540M/611M [00:06<00:00, 82.3MB/s]
 90%|███████████████████████████████████    | 549M/611M [00:06<00:00, 84.0MB/s]
 91%|███████████████████████████████████▌   | 558M/611M [00:06<00:00, 85.1MB/s]
 93%|████████████████████████████████████▏  | 566M/611M [00:07<00:00, 85.9MB/s]
 94%|████████████████████████████████████▋  | 575M/611M [00:07<00:00, 86.3MB/s]
 96%|█████████████████████████████████████▎ | 584M/611M [00:07<00:00, 86.3MB/s]
 97%|█████████████████████████████████████▊ | 593M/611M [00:07<00:00, 87.0MB/s]
 98%|██████████████████████████████████████▍| 601M/611M [00:07<00:00, 87.2MB/s]
100%|██████████████████████████████████████▉| 610M/611M [00:07<00:00, 87.3MB/s]
  0%|                                               | 0.00/611M [00:00<?, ?B/s]
100%|███████████████████████████████████████| 611M/611M [00:00<00:00, 2.36TB/s]
Download complete in 18s (582.2 MB)
Opening raw data file /home/circleci/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
    Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
Ready.
111 events found on stim channel STI 014
Event IDs: [1]
Reading 0 ... 269399  =      0.000 ...   897.077 secs...
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 35 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 8.75 Hz (-6 dB cutoff frequency: 39.38 Hz)
- Filter length: 993 samples (3.307 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.6s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    1.4s
Not setting metadata
111 matching events found
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 111 events and 251 original time points ...
0 bad epochs dropped

Assessing connectivity in non-evoked data

We will first demonstrate how connectivity can be assessed from non-evoked data. In this example, we use data from the pre-trial period of [-1.5, -0.5] seconds. We compute Fourier coefficients of the data using the compute_psd() method with output="complex" (note that this requires mne >= 1.8).

Next, we pass these coefficients to spectral_connectivity_epochs() to compute connectivity using the imaginary part of coherency (imcoh). Our indices specify that connectivity should be computed between all pairs of channels.

# Compute Fourier coefficients for pre-trial data
fmin, fmax = 3, 23
pretrial_coeffs = epochs.compute_psd(
    fmin=fmin, fmax=fmax, tmin=None, tmax=-0.5, output="complex"
)
freqs = pretrial_coeffs.freqs

# Compute connectivity for pre-trial data
indices = np.tril_indices(epochs.info["nchan"], k=-1)  # all-to-all connectivity
pretrial_con = spectral_connectivity_epochs(
    pretrial_coeffs, method="imcoh", indices=indices
)

    Using multitaper spectrum estimation with 7 DPSS windows
Connectivity computation...
    frequencies: 4.0Hz..22.8Hz (20 points)
    computing connectivity for 20706 connections
    the following metrics will be computed: Imaginary Coherence
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
[Connectivity computation done]

Next, we generate the surrogate data by passing the Fourier coefficients into the make_surrogate_data() function. To get a reliable estimate of the baseline connectivity, we perform this shuffling procedure \(\text{n}_{\text{shuffle}}\) times, producing \(\text{n}_{\text{shuffle}}\) surrogate datasets. We can then iterate over these shuffles and compute the connectivity for each one.

# Generate surrogate data
n_shuffles = 100  # recommended is >= 1,000; limited here to reduce compute time
pretrial_surrogates = make_surrogate_data(
    pretrial_coeffs, n_shuffles=n_shuffles, rng_seed=44
)

# Compute connectivity for surrogate data
surrogate_con = []
for shuffle_i, surrogate in enumerate(pretrial_surrogates, 1):
    print(f"Computing connectivity for shuffle {shuffle_i} of {n_shuffles}")
    surrogate_con.append(
        spectral_connectivity_epochs(
            surrogate, method="imcoh", indices=indices, verbose=False
        )
    )

Computing connectivity for shuffle 1 of 100
Computing connectivity for shuffle 2 of 100
Computing connectivity for shuffle 3 of 100
Computing connectivity for shuffle 4 of 100
Computing connectivity for shuffle 5 of 100
Computing connectivity for shuffle 6 of 100
Computing connectivity for shuffle 7 of 100
Computing connectivity for shuffle 8 of 100
Computing connectivity for shuffle 9 of 100
Computing connectivity for shuffle 10 of 100
Computing connectivity for shuffle 11 of 100
Computing connectivity for shuffle 12 of 100
Computing connectivity for shuffle 13 of 100
Computing connectivity for shuffle 14 of 100
Computing connectivity for shuffle 15 of 100
Computing connectivity for shuffle 16 of 100
Computing connectivity for shuffle 17 of 100
Computing connectivity for shuffle 18 of 100
Computing connectivity for shuffle 19 of 100
Computing connectivity for shuffle 20 of 100
Computing connectivity for shuffle 21 of 100
Computing connectivity for shuffle 22 of 100
Computing connectivity for shuffle 23 of 100
Computing connectivity for shuffle 24 of 100
Computing connectivity for shuffle 25 of 100
Computing connectivity for shuffle 26 of 100
Computing connectivity for shuffle 27 of 100
Computing connectivity for shuffle 28 of 100
Computing connectivity for shuffle 29 of 100
Computing connectivity for shuffle 30 of 100
Computing connectivity for shuffle 31 of 100
Computing connectivity for shuffle 32 of 100
Computing connectivity for shuffle 33 of 100
Computing connectivity for shuffle 34 of 100
Computing connectivity for shuffle 35 of 100
Computing connectivity for shuffle 36 of 100
Computing connectivity for shuffle 37 of 100
Computing connectivity for shuffle 38 of 100
Computing connectivity for shuffle 39 of 100
Computing connectivity for shuffle 40 of 100
Computing connectivity for shuffle 41 of 100
Computing connectivity for shuffle 42 of 100
Computing connectivity for shuffle 43 of 100
Computing connectivity for shuffle 44 of 100
Computing connectivity for shuffle 45 of 100
Computing connectivity for shuffle 46 of 100
Computing connectivity for shuffle 47 of 100
Computing connectivity for shuffle 48 of 100
Computing connectivity for shuffle 49 of 100
Computing connectivity for shuffle 50 of 100
Computing connectivity for shuffle 51 of 100
Computing connectivity for shuffle 52 of 100
Computing connectivity for shuffle 53 of 100
Computing connectivity for shuffle 54 of 100
Computing connectivity for shuffle 55 of 100
Computing connectivity for shuffle 56 of 100
Computing connectivity for shuffle 57 of 100
Computing connectivity for shuffle 58 of 100
Computing connectivity for shuffle 59 of 100
Computing connectivity for shuffle 60 of 100
Computing connectivity for shuffle 61 of 100
Computing connectivity for shuffle 62 of 100
Computing connectivity for shuffle 63 of 100
Computing connectivity for shuffle 64 of 100
Computing connectivity for shuffle 65 of 100
Computing connectivity for shuffle 66 of 100
Computing connectivity for shuffle 67 of 100
Computing connectivity for shuffle 68 of 100
Computing connectivity for shuffle 69 of 100
Computing connectivity for shuffle 70 of 100
Computing connectivity for shuffle 71 of 100
Computing connectivity for shuffle 72 of 100
Computing connectivity for shuffle 73 of 100
Computing connectivity for shuffle 74 of 100
Computing connectivity for shuffle 75 of 100
Computing connectivity for shuffle 76 of 100
Computing connectivity for shuffle 77 of 100
Computing connectivity for shuffle 78 of 100
Computing connectivity for shuffle 79 of 100
Computing connectivity for shuffle 80 of 100
Computing connectivity for shuffle 81 of 100
Computing connectivity for shuffle 82 of 100
Computing connectivity for shuffle 83 of 100
Computing connectivity for shuffle 84 of 100
Computing connectivity for shuffle 85 of 100
Computing connectivity for shuffle 86 of 100
Computing connectivity for shuffle 87 of 100
Computing connectivity for shuffle 88 of 100
Computing connectivity for shuffle 89 of 100
Computing connectivity for shuffle 90 of 100
Computing connectivity for shuffle 91 of 100
Computing connectivity for shuffle 92 of 100
Computing connectivity for shuffle 93 of 100
Computing connectivity for shuffle 94 of 100
Computing connectivity for shuffle 95 of 100
Computing connectivity for shuffle 96 of 100
Computing connectivity for shuffle 97 of 100
Computing connectivity for shuffle 98 of 100
Computing connectivity for shuffle 99 of 100
Computing connectivity for shuffle 100 of 100

We can plot the all-to-all connectivity of the pre-trial data against the surrogate data, averaged over all shuffles. This shows a strong degree of coupling in the alpha band (~8-12 Hz), with weaker coupling in the lower range of the beta band (~13-20 Hz). A simple visual inspection shows that connectivity in the alpha and beta bands are above the baseline level of connectivity estimated from the surrogate data. However, we need to confirm this statistically.

# Plot pre-trial vs. surrogate connectivity
fig, ax = plt.subplots(1, 1)
ax.plot(
    freqs,
    np.abs([surrogate.get_data() for surrogate in surrogate_con]).mean(axis=(0, 1)),
    linestyle="--",
    label="Surrogate",
)
ax.plot(freqs, np.abs(pretrial_con.get_data()).mean(axis=0), label="Original")
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Connectivity (A.U.)")
ax.set_title("All-to-all connectivity | Pre-trial ")
ax.legend()

<matplotlib.legend.Legend object at 0x7f590ea57ac0>

Assessing the statistical significance of our connectivity estimates can be done with the following simple procedure [2]

\(p=\LARGE{\frac{\Sigma_{s=1}^Sc_s}{S}}\) ,

\(c_s=\{1\text{ if }\text{Con}\leq\text{Con}_{\text{s}}\text{ },\text{ }0 \text{ if otherwise }\) ,

where: \(p\) is our p-value; \(s\) is a given shuffle iteration of \(S\) total shuffles; and \(c\) is a binary indicator of whether the true connectivity, \(\text{Con}\), is greater than the surrogate connectivity, \(\text{Con}_{\text{s}}\), for a given shuffle.

Note that for connectivity methods which produce negative scores (e.g., imaginary part of coherency, time-reversed Granger causality, etc…), you should take the absolute values before testing. Similar adjustments should be made for methods that produce scores centred around non-zero values (e.g., 0.5 for directed phase lag index).

Below, we determine the statistical significance of connectivity in the lower beta band. We simplify this by averaging over all connections and corresponding frequency bins. We could of course also test the significance of each connection, each frequency bin, or other frequency bands such as the alpha band. Naturally, any tests involving multiple connections, frequencies, and/or times should be corrected for multiple comparisons.

The test confirms our visual inspection, showing that connectivity in the lower beta band is significantly above the baseline level of connectivity at an alpha of 0.05, which we can take as evidence of genuine interactions in this frequency band.

# Find indices of lower beta frequencies
beta_freqs = np.where((freqs >= 13) & (freqs <= 20))[0]

# Compute lower beta connectivity for pre-trial data (average connections and freqs)
beta_con_pretrial = np.abs(pretrial_con.get_data()[:, beta_freqs]).mean(axis=(0, 1))

# Compute lower beta connectivity for surrogate data (average connections and freqs)
beta_con_surrogate = np.abs(
    [surrogate.get_data()[:, beta_freqs] for surrogate in surrogate_con]
).mean(axis=(1, 2))

# Compute p-value for pre-trial lower beta coupling
p_val = np.sum(beta_con_pretrial <= beta_con_surrogate) / n_shuffles
print(f"P = {p_val:.2f}")

P = 0.00

Assessing connectivity in evoked data

When generating surrogate data, it is important to distinguish non-evoked data (e.g., resting-state, pre/inter-trial data) from evoked data (where a stimulus is presented or an action performed at a set time during each epoch). Critically, evoked data contains a temporal structure that is consistent across epochs, and thus shuffling epochs across channels will fail to adequately disrupt the covariance structure.

Any connectivity estimates will therefore overestimate the baseline connectivity in your data, increasing the likelihood of type II errors (see the Notes section of make_surrogate_data() for more information, and see the section Generating surrogate connectivity from inappropriate data for a demonstration).

In cases where you want to assess connectivity in evoked data, you can use surrogates generated from non-evoked data (of the same subject). Here we do just that, comparing connectivity estimates from the pre-trial surrogates to the evoked, post-stimulus response ([0, 1] second).

Again, there is pronounced alpha coupling (stronger than in the pre-trial data) and weaker beta coupling, both of which appear to be above the baseline level of connectivity.

# Compute Fourier coefficients for post-stimulus data
poststim_coeffs = epochs.compute_psd(
    fmin=fmin, fmax=fmax, tmin=0, tmax=None, output="complex"
)

# Compute connectivity for post-stimulus data
poststim_con = spectral_connectivity_epochs(
    poststim_coeffs, method="imcoh", indices=indices
)

# Plot post-stimulus vs. (pre-trial) surrogate connectivity
fig, ax = plt.subplots(1, 1)
ax.plot(
    freqs,
    np.abs([surrogate.get_data() for surrogate in surrogate_con]).mean(axis=(0, 1)),
    linestyle="--",
    label="Surrogate",
)
ax.plot(freqs, np.abs(poststim_con.get_data()).mean(axis=0), label="Original")
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Connectivity (A.U.)")
ax.set_title("All-to-all connectivity | Post-stimulus")
ax.legend()

    Using multitaper spectrum estimation with 7 DPSS windows
Connectivity computation...
    frequencies: 4.0Hz..22.8Hz (20 points)
    computing connectivity for 20706 connections
    the following metrics will be computed: Imaginary Coherence
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
[Connectivity computation done]

<matplotlib.legend.Legend object at 0x7f58a89d7b20>

This is also confirmed by statistical testing, with connectivity in the lower beta band being significantly above the baseline level of connectivity. Thus, using surrogate connectivity estimates from non-evoked data provides a reliable baseline for assessing connectivity in evoked data.

# Compute lower beta connectivity for post-stimulus data (average connections and freqs)
beta_con_poststim = np.abs(poststim_con.get_data()[:, beta_freqs]).mean(axis=(0, 1))

# Compute p-value for post-stimulus lower beta coupling
p_val = np.sum(beta_con_poststim <= beta_con_surrogate) / n_shuffles
print(f"P = {p_val:.2f}")

P = 0.00

Generating surrogate connectivity from inappropriate data

We discussed above how surrogates generated from evoked data risk overestimating the degree of baseline connectivity. We demonstrate this below by generating surrogates from the post-stimulus data.

# Generate surrogates from evoked data
poststim_surrogates = make_surrogate_data(
    poststim_coeffs, n_shuffles=n_shuffles, rng_seed=44
)

# Compute connectivity for evoked surrogate data
bad_surrogate_con = []
for shuffle_i, surrogate in enumerate(poststim_surrogates, 1):
    print(f"Computing connectivity for shuffle {shuffle_i} of {n_shuffles}")
    bad_surrogate_con.append(
        spectral_connectivity_epochs(
            surrogate, method="imcoh", indices=indices, verbose=False
        )
    )

Computing connectivity for shuffle 1 of 100
Computing connectivity for shuffle 2 of 100
Computing connectivity for shuffle 3 of 100
Computing connectivity for shuffle 4 of 100
Computing connectivity for shuffle 5 of 100
Computing connectivity for shuffle 6 of 100
Computing connectivity for shuffle 7 of 100
Computing connectivity for shuffle 8 of 100
Computing connectivity for shuffle 9 of 100
Computing connectivity for shuffle 10 of 100
Computing connectivity for shuffle 11 of 100
Computing connectivity for shuffle 12 of 100
Computing connectivity for shuffle 13 of 100
Computing connectivity for shuffle 14 of 100
Computing connectivity for shuffle 15 of 100
Computing connectivity for shuffle 16 of 100
Computing connectivity for shuffle 17 of 100
Computing connectivity for shuffle 18 of 100
Computing connectivity for shuffle 19 of 100
Computing connectivity for shuffle 20 of 100
Computing connectivity for shuffle 21 of 100
Computing connectivity for shuffle 22 of 100
Computing connectivity for shuffle 23 of 100
Computing connectivity for shuffle 24 of 100
Computing connectivity for shuffle 25 of 100
Computing connectivity for shuffle 26 of 100
Computing connectivity for shuffle 27 of 100
Computing connectivity for shuffle 28 of 100
Computing connectivity for shuffle 29 of 100
Computing connectivity for shuffle 30 of 100
Computing connectivity for shuffle 31 of 100
Computing connectivity for shuffle 32 of 100
Computing connectivity for shuffle 33 of 100
Computing connectivity for shuffle 34 of 100
Computing connectivity for shuffle 35 of 100
Computing connectivity for shuffle 36 of 100
Computing connectivity for shuffle 37 of 100
Computing connectivity for shuffle 38 of 100
Computing connectivity for shuffle 39 of 100
Computing connectivity for shuffle 40 of 100
Computing connectivity for shuffle 41 of 100
Computing connectivity for shuffle 42 of 100
Computing connectivity for shuffle 43 of 100
Computing connectivity for shuffle 44 of 100
Computing connectivity for shuffle 45 of 100
Computing connectivity for shuffle 46 of 100
Computing connectivity for shuffle 47 of 100
Computing connectivity for shuffle 48 of 100
Computing connectivity for shuffle 49 of 100
Computing connectivity for shuffle 50 of 100
Computing connectivity for shuffle 51 of 100
Computing connectivity for shuffle 52 of 100
Computing connectivity for shuffle 53 of 100
Computing connectivity for shuffle 54 of 100
Computing connectivity for shuffle 55 of 100
Computing connectivity for shuffle 56 of 100
Computing connectivity for shuffle 57 of 100
Computing connectivity for shuffle 58 of 100
Computing connectivity for shuffle 59 of 100
Computing connectivity for shuffle 60 of 100
Computing connectivity for shuffle 61 of 100
Computing connectivity for shuffle 62 of 100
Computing connectivity for shuffle 63 of 100
Computing connectivity for shuffle 64 of 100
Computing connectivity for shuffle 65 of 100
Computing connectivity for shuffle 66 of 100
Computing connectivity for shuffle 67 of 100
Computing connectivity for shuffle 68 of 100
Computing connectivity for shuffle 69 of 100
Computing connectivity for shuffle 70 of 100
Computing connectivity for shuffle 71 of 100
Computing connectivity for shuffle 72 of 100
Computing connectivity for shuffle 73 of 100
Computing connectivity for shuffle 74 of 100
Computing connectivity for shuffle 75 of 100
Computing connectivity for shuffle 76 of 100
Computing connectivity for shuffle 77 of 100
Computing connectivity for shuffle 78 of 100
Computing connectivity for shuffle 79 of 100
Computing connectivity for shuffle 80 of 100
Computing connectivity for shuffle 81 of 100
Computing connectivity for shuffle 82 of 100
Computing connectivity for shuffle 83 of 100
Computing connectivity for shuffle 84 of 100
Computing connectivity for shuffle 85 of 100
Computing connectivity for shuffle 86 of 100
Computing connectivity for shuffle 87 of 100
Computing connectivity for shuffle 88 of 100
Computing connectivity for shuffle 89 of 100
Computing connectivity for shuffle 90 of 100
Computing connectivity for shuffle 91 of 100
Computing connectivity for shuffle 92 of 100
Computing connectivity for shuffle 93 of 100
Computing connectivity for shuffle 94 of 100
Computing connectivity for shuffle 95 of 100
Computing connectivity for shuffle 96 of 100
Computing connectivity for shuffle 97 of 100
Computing connectivity for shuffle 98 of 100
Computing connectivity for shuffle 99 of 100
Computing connectivity for shuffle 100 of 100

Plotting the post-stimulus connectivity against the estimates from the non-evoked and evoked surrogate data, we see that the evoked surrogate data greatly overestimates the baseline connectivity in the alpha band.

Although in this case the alpha connectivity was still far above the baseline from the evoked surrogates, this will not always be the case, and you can see how this risks false negative assessments that connectivity is not significantly different from baseline.

# Plot post-stimulus vs. evoked and non-evoked surrogate connectivity
fig, ax = plt.subplots(1, 1)
ax.plot(
    freqs,
    np.abs([surrogate.get_data() for surrogate in surrogate_con]).mean(axis=(0, 1)),
    linestyle="--",
    label="Surrogate (pre-stimulus)",
)
ax.plot(
    freqs,
    np.abs([surrogate.get_data() for surrogate in bad_surrogate_con]).mean(axis=(0, 1)),
    color="C3",
    linestyle="--",
    label="Surrogate (post-stimulus)",
)
ax.plot(
    freqs, np.abs(poststim_con.get_data()).mean(axis=0), color="C1", label="Original"
)
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Connectivity (A.U.)")
ax.set_title("All-to-all connectivity | Post-stimulus")
ax.legend()

<matplotlib.legend.Legend object at 0x7f590a0e32b0>

Assessing connectivity on a group-level

While our focus here has been on assessing the significance of connectivity on a single recording-level, we may also want to determine whether group-level connectivity estimates are significantly different from baseline. For this, we can generate surrogates and estimate connectivity alongside the original signals for each piece of data.

There are multiple ways to assess the statistical significance. For example, we can compute p-values for each piece of data using the approach above and combine them for the nested data (e.g., across recordings, subjects, etc…) using Stouffer’s method [3].

Alternatively, we could take the average of the surrogate connectivity estimates across all shuffles for each piece of data and compare them to the original connectivity estimates in a paired test. The scipy.stats and mne.stats modules have many such tools for testing this, e.g., scipy.stats.ttest_1samp(), mne.stats.permutation_t_test(), etc…

Altogether, surrogate connectivity estimates are a powerful tool for assessing the significance of connectivity estimates, both on a single recording- and group-level.

References

[1]

Martin Vinck, Marijn van Wingerden, Thilo Womelsdorf, Pascal Fries, and Cyriel M.A. Pennartz. The pairwise phase consistency: a bias-free measure of rhythmic neuronal synchronization. NeuroImage, 51(1):112–122, 2010. doi:10.1016/j.neuroimage.2010.01.073.

[2] (1,2)

Franziska Pellegrini, Arnaud Delorme, Vadim Nikulin, and Stefan Haufe. Identifying good practices for detecting inter-regional linear functional connectivity from EEG. NeuroImage, 277:120218, 2023. doi:10.1016/j.neuroimage.2023.120218.

[3]

Irene Dowding and Stefan Haufe. Powerful statistical inference for nested data using sufficient summary statistics. Frontiers in Human Neuroscience, 12:103, 2018. doi:10.3389/fnhum.2018.00103.