Example demonstrating difference between spectral connectivity over Epochs or over time

[adam2392]

With the inclusion of #67 spectral_connectivity_epochs and spectral_connectivity_time are both available to estimate frequency-dependent connectivity.

It would be helpful to augment the existing documentation with an example demonstrating the differences between the two.

Anyone in the community willing to help out with this? Would definitely be considered as a co-author in a JOSS/Methods publication for mne-connectivity.

[Div12345]

Hey @adam2392 , I could make an attempt at this.

[adam2392]

Sure go for it! I think ideally demonstrating this with a simulation + real dataset would be valuable.

[Avoide]

@adam2392 In #17 I also tested the difference between spectral connectivity over trials/Epochs vs time. I have adapted the code to use spectral_connectivity_epochs vs spectral_connectivity_time. It is quite a simple simulation, the first data option uses the same repeated trial (hence spectral connectivity over Epochs is 1, but not over time since it is just random). The second data generation option just creates 10Hz sinus waves with different random phases, however they are created differently for each epoch. This means the connectivity over time is 1, but not for con over epochs.

It might be useful for inspiration.

# Libraries
import numpy as np
import mne
from mne_connectivity import spectral_connectivity_time
from mne_connectivity import spectral_connectivity_epochs

# Generate some data
n_epochs = 5
n_channels = 3
n_times = 1000
sfreq = 250

# Choose what kind of data should be used for evaluation
# 0 = same trial repeated in all epochs - we expect high connectivity over Epochs
# 1 = 10Hz sinus waves with random phase differences in each channel and epoch - we expect high connectivity over time
data_option = 0
np.random.seed(42)
data = np.random.rand(n_epochs, n_channels, n_times)

if data_option == 0:
    print("Data option 0: Expect high Con over Epochs")
    # Make all 5 epochs the same trial to show difference between connectivity
    # over time and over trials. Here we expect Con over trials = 1
    for i in range(n_epochs):
        data[i] = data[0]
elif data_option == 1:
    print("Data option 0: Expect high 10Hz Con over Time")
    # Generate 10Hz sinus waves to show difference between connectivity
    # over time and over trials. Here we expect Con over time = 1
    for i in range(n_epochs):
        for c in range(n_channels):
            wave_freq = 10
            epoch_len = n_times/sfreq
            # Introduce random phase for each channel
            phase = np.random.rand(1)*10
            # Generate sinus wave
            x = np.linspace(-wave_freq*epoch_len*np.pi+phase,
                            wave_freq*epoch_len*np.pi+phase,n_times)
            data[i,c] = np.squeeze(np.sin(x))
else:
    print("Data_option value chosen is invalid")

# Create epochs object for compatibility
ch_names = ["T1","T2","T3"] # random names
info = mne.create_info(ch_names, sfreq, ch_types="eeg")
data_epoch = mne.EpochsArray(data,info)

# Visualize the data
data_epoch.plot(scalings=1)

# Calculate coh over epochs/trials
con_Epoch = spectral_connectivity_epochs(data_epoch, method="coh",
    mode="cwt_morlet", sfreq=sfreq, cwt_freqs=np.array([10]))

con_Epoch = con_Epoch.get_data(output="dense")
# Avg over time
con_Epoch = np.mean(con_Epoch,axis=-1)
print("10Hz Coh over Epochs")
print(con_Epoch[:,:,0]) # all 1, since it is the same trial

# Calculate coh over time
con_time = spectral_connectivity_time(data_epoch, method="coh",
    mode="cwt_morlet", sfreq=sfreq, freqs=10)

con_time = con_time.get_data()
# Avg over time
con_time = np.mean(con_time,axis=-1)
# Avg over epochs
con_time = np.mean(con_time,axis=0)
print("10Hz Coh over time")
print(con_time[:,:,0]) # all 1, since it is the same trial

# Using data_option 1
data_option = 1
np.random.seed(42)
data = np.random.rand(n_epochs, n_channels, n_times)

if data_option == 0:
    print("Data option 0: Expect high Con over Epochs")
    # Make all 5 epochs the same trial to show difference between connectivity
    # over time and over trials. Here we expect Con over trials = 1
    for i in range(n_epochs):
        data[i] = data[0]
elif data_option == 1:
    print("Data option 0: Expect high 10Hz Con over Time")
    # Generate 10Hz sinus waves to show difference between connectivity
    # over time and over trials. Here we expect Con over time = 1
    for i in range(n_epochs):
        for c in range(n_channels):
            wave_freq = 10
            epoch_len = n_times/sfreq
            # Introduce random phase for each channel
            phase = np.random.rand(1)*10
            # Generate sinus wave
            x = np.linspace(-wave_freq*epoch_len*np.pi+phase,
                            wave_freq*epoch_len*np.pi+phase,n_times)
            data[i,c] = np.squeeze(np.sin(x))
else:
    print("Data_option value chosen is invalid")

# Create epochs object for compatibility
ch_names = ["T1","T2","T3"] # random names
info = mne.create_info(ch_names, sfreq, ch_types="eeg")
data_epoch = mne.EpochsArray(data,info)

# Visualize the data
data_epoch.plot(scalings=2)

# Calculate coh over epochs/trials
con_Epoch = spectral_connectivity_epochs(data_epoch, method="coh",
    mode="cwt_morlet", sfreq=sfreq, cwt_freqs=np.array([10]))

con_Epoch = con_Epoch.get_data(output="dense")
# Avg over time
con_Epoch = np.mean(con_Epoch,axis=-1)
print("10Hz Coh over Epochs")
print(con_Epoch[:,:,0]) # all 1, since it is the same trial

# Calculate coh over time
con_time = spectral_connectivity_time(data_epoch, method="coh",
    mode="cwt_morlet", sfreq=sfreq, freqs=10)

con_time = con_time.get_data()
# Avg over time
con_time = np.mean(con_time,axis=-1)
# Avg over epochs
con_time = np.mean(con_time,axis=0)
print("10Hz Coh over time")
print(con_time[:,:,0]) # all 1, since it is the same trial

[adam2392]

Hi, it seems that the issues raised in the discussion above have been resolved now. However, just for transparency, we are still interested in accepting a contribution regarding the example docs. Thanks for all the discussion so far!

[Avoide]

If there are no one actively working on an example doc for demonstrating the difference of connectivity over time or over epochs, then I could try and make one. The simulation should be straightforward, but I'm a bit in doubt about what real world dataset to use, so if anyone have some ideas it would be helpful.

[adam2392]

Perhaps any on openneuro, if we can easily clip the dataset, or perhaps one from the existing MNE-Python, or MNE-connectivity tutorials would be good candidates perhaps :)

[xavi0488]

@adam2392 : Thank you so much for your great work MNE. @Avoide : Thanks for your simulation code, it helps me a lot. I modified it as follow because I found some errors and incorrect. Here I used PLI instead of Coh.

import numpy as np
import mne
from mne_connectivity import spectral_connectivity_time
from mne_connectivity import spectral_connectivity_epochs

# Generate some data
n_epochs = 5
n_channels = 3
n_times = 1000
sfreq = 250

# CASE 1-----------------------------------------------------------------------
np.random.seed(42)
data = np.random.rand(n_epochs, n_channels, n_times)

for i in range(n_epochs):
    data[i] = data[0]

ch_names = ["T1","T2","T3"] # random names
info = mne.create_info(ch_names, sfreq, ch_types="eeg")
data_epoch = mne.EpochsArray(data,info)
data_epoch.plot(scalings=1) # Visualize the data

# Using spectral_connectivity_epochs: phase difference is quantified over epochs at each specific time point of epoch
con_Epoch = spectral_connectivity_epochs(data_epoch, method="pli",
    mode="cwt_morlet", sfreq=sfreq, cwt_freqs=np.array([10]))
con_Epoch = con_Epoch.get_data(output="dense")
con_Epoch = np.mean(con_Epoch[:,:,0,:],axis=-1) # avg over times

# Using spectral_connectivity_time: phase difference is quantified over all time points of each epoch
con_time = spectral_connectivity_time(data_epoch, method="pli",
    mode="cwt_morlet", sfreq=sfreq, freqs=10)
con_time = con_time.get_data()
con_time = np.mean(con_time[:,:,0],axis=0) # avg over epochs

# CASE 2-----------------------------------------------------------------------
np.random.seed(42)
data = np.random.rand(n_epochs, n_channels, n_times)

for i in range(n_epochs):
    for c in range(n_channels):
        wave_freq = 10
        epoch_len = n_times/sfreq
        # Introduce random phase for each channel
        phase = np.random.rand(1)*10
        # Generate sinus wave
        x = np.linspace(-wave_freq*epoch_len*np.pi+phase,
                            wave_freq*epoch_len*np.pi+phase,n_times)
        data[i,c] = np.squeeze(np.sin(x))

ch_names = ["T1","T2","T3"] # random names
info = mne.create_info(ch_names, sfreq, ch_types="eeg")
data_epoch = mne.EpochsArray(data,info)
data_epoch.plot(scalings=2) # Visualize the data

# Using spectral_connectivity_epochs: phase difference is quantified over epochs at each specific time point of epoch
con_Epoch = spectral_connectivity_epochs(data_epoch, method="pli",
    mode="cwt_morlet", sfreq=sfreq, cwt_freqs=np.array([10]))
con_Epoch = con_Epoch.get_data(output="dense")
con_Epoch = np.mean(con_Epoch[:,:,0,:],axis=-1) # avg over times

# Using spectral_connectivity_time: phase difference is quantified over all time points of each epoch
con_time = spectral_connectivity_time(data_epoch, method="pli",
    mode="cwt_morlet", sfreq=sfreq, freqs=10)
con_time = con_time.get_data()
con_time = np.mean(con_time[:,:,0],axis=0) # avg over epochs

[adam2392]

@adam2392 : Thank you so much for your great work MNE. @Avoide : Thanks for your simulation code, it helps me a lot. I modified it as follow because I found some errors and incorrect. Here I used PLI instead of Coh.

What are the errors/issue? Do you mind helping to summarizing it here? And what version of mne are you using? E.g. the output of mne sys_info

[xavi0488]

Dear @adam2392 , Your tool is great, it has not any errors. I use MNE 1.0.3. The errors I found came from the code of @Avoide when he prints the result. Moreover, some results should not 1 but he remarked it is 1. This may confuse some readers.

[Avoide]

@TrinhThanhTung ahh yes I see what you mean with the comments that stated it was 1 when it shouldn't be. Thanks for the heads-up. I copied the code and just changed the "Data setting" variable, but forgot to change the comments in my previous post. The term "case" in the code you posted is probably more intuitive. I'll keep it in mind when writing the example doc.

[xavi0488]

@Avoide A lot of thanks to your codes!!!

[xavi0488]

Dear @Avoide, @adam2392 , Could you please help to answer this question:

I have 120 seconds of resting-state EEG (sampling rate=500 Hz). I segment it into 40 non-overlap 3 seconds, making the data of the shape (40,30,1500). Here 30 is the number of EEG electrodes. I think in this case "spectral_connectivity_time" is the right selection to calculate the connectivity of electrode pairs. If I want to calculate the alpha band (8~13 Hz) connectivity, how to set the parameters": freqs, fmin, fmax in your "spectral_connectivity_time" function?

Thank you so much!

[Avoide]

@TrinhThanhTung Yes for resting-state EEG, I would choose spectral_connectivity_time.

You could try:

Freq_Bands = {"alpha": [8.0, 13.0]}
min_freq = np.min(list(Freq_Bands.values()))
max_freq = np.max(list(Freq_Bands.values()))

freqs = np.linspace(min_freq, max_freq, int((max_freq - min_freq) * 4 + 1))
fmin = tuple([list(Freq_Bands.values())[f][0] for f in range(len(Freq_Bands))])
fmax = tuple([list(Freq_Bands.values())[f][1] for f in range(len(Freq_Bands))])

This code allows you to change the freq band in the dictionary, but also add more frequency bands, e.g. theta or beta.

Freq_Bands = {"theta": [4.0, 8.0],
              "alpha": [8.0, 13.0],
              "beta": [13.0, 30.0]}

And the freqs, fmin and fmax will change accordingly.

[xavi0488]

@Avoide I really appreciate the time you put into helping me with a very clear answer. I tried and the code works well for all bands, except the delta band (1~4 Hz).

It raises this error: *** ValueError: At least one value in n_cycles corresponds to a wavelet longer than the signal. Use less cycles, higher frequencies, or longer epochs.

I used: method='pli', mode='cwt_morlet'.

Actually, I wrote my own code (quite similar to your code) and met above error when I try with delta band first. Can you suggest what problem may be?

[adam2392]

ValueError: At least one value in n_cycles corresponds to a wavelet longer than the signal. Use less cycles, higher frequencies, or longer epochs.

This most likely means your data is not long enough.

Btw I think pli is implemented inside mne-connectivity, so you don't need to write your own implementation anymore?

Btw @TrinhThanhTung we have a great Discourse webpage for answering usage questions exactly like this: https://mne.discourse.group. Perhaps you'll find more fruitful discussion there since it's visible to everyone + searchable (so maybe someone previously had this issue).

[Avoide]

@TrinhThanhTung The error is due to the short duration of the epochs. If you have 3s epochs, but want to investigate 1 Hz freq, then you can only use 3 cycles. Otherwise the 1 Hz wavelet will be longer than the signal. It pertains to the trade-off between time resolution and frequency resolution. The better freq resolution you want, the worse time resolution.

[xavi0488]

@adam2392 Thank you so much for your recommendation. @Avoide Oh, I see. You had a clear explain. Thanks a lot!

[ruuskas]

@TrinhThanhTung I recommend checking out the docs for mne.time_frequency.tfr_array_morlet for better understanding of the limitations of the time-frequency analysis underlying the connectivity computation.

[xavi0488]

@ruuskas Thank you so much for your recommendation.