Time-Resolved Spectral Connectivity PLI implementation

[christian-oreilly]

Describe the bug

I think the implementation of the PLI metric is wrong. The standard definition is:

PLI=|mean(sign(difference_in_instantaneous_phase))|

Note that for a signal of N time samples, this PLI value is obtained as the average of N values. The MNE implementation of this metric follows the implementation of Ortiz and instead of computing this metric in the time domain (e.g. using the Hilbert transform to estimate the instantaneous phase) computes it in the frequency domain using the cross-spectrum. An issue arise when this passage to the frequency domain is used to analyze the PLI results at the individual frequency bands. For example, the default behavior of spectral_connectivity() is faverage=False, which therefore report a value for each frequency. In this case, the PLI estimate that was the absolute-mean-value of N samples become the absolute-mean-value of 1 sample... so the PLI estimates reduce to 0 or 1 (where it was originally a value taken from the {0/N, 1/N, 2/N, ..., N/N}). Now, if we might think that we would recover the correct values using faverage=True, but we don't because this computes the mean-absolute-value, not the absolute-mean-value... Here is a part of the spectral_connectivity() code:

    for method, n_args in zip(con_methods, n_comp_args):
        # future estimators will need to be handled here
        if n_args == 3:
            # compute all scores at once
            method.compute_con(slice(0, n_cons), n_epochs)
        elif n_args == 5:
            # compute scores block-wise to save memory
            for i in range(0, n_cons, block_size):
                con_idx = slice(i, i + block_size)
                psd_xx = psd[idx_map[0][con_idx]]
                psd_yy = psd[idx_map[1][con_idx]]
                method.compute_con(con_idx, n_epochs, psd_xx, psd_yy)
        else:
            raise RuntimeError('This should never happen.')

        # get the connectivity scores
        this_con = method.con_scores

        if this_con.shape[0] != n_cons:
            raise ValueError('First dimension of connectivity scores must be '
                             'the same as the number of connections')
        if faverage:
            if this_con.shape[1] != n_freqs:
                raise ValueError('2nd dimension of connectivity scores must '
                                 'be the same as the number of frequencies')
            con_shape = (n_cons, n_bands) + this_con.shape[2:]
            this_con_bands = np.empty(con_shape, dtype=this_con.dtype)
            for band_idx in range(n_bands):
                this_con_bands[:, band_idx] =\
                    np.mean(this_con[:, freq_idx_bands[band_idx]], axis=1)
            this_con = this_con_bands

Note that the faverage part comes after the call to compute_con() which, in the case of _PLIEst already applied the absolute operator. So, whereas in the definition, you average N values taken in {-1, 0, 1}, here you average N values taken in {0, 1}.

If I am not mistaken, the function "appears" to work when using epochs because then it does the mean(signs(X)) across the M epochs, before computing the absolute value. But remember that the stuff in X was initially N values saying whether the time-lag at every N time samples was positive or negative. Now X values are still saying whether the phase-difference is positive or negative, but you traded the ability to check it at every time point to have it localized in frequencies (uncertainty principle), so you have just one value per epochs and you average it across M epochs. I might have interesting properties on its own, but I don't think it is still fair to call this a PLI measure, i.e. average_across_epochs(abs(average_across_x)) != average_across_x(abs(average_across_epochs)) (where x can be either time or frequency, depending on whether you transform from temporal representation to frequency).

[christian-oreilly]

As an afterthought (why do we need to press "submit" to have these?), I guess that when resting-state continuous signals are "chunked" in epochs, the "epoch dimension" has a "temporal definition" to it so it still makes some sense in that case. But I don't think that it makes sense for epochs that have a meaning of "repetition" like in ERP studies... In any case, using the standard PLI definition separately on every epoch and then averaging these PLI values across epochs would give a different value.

Also, it is definitely problematic when someone wants to compute the standard PLI value between two-signals (no epoching).

Anyway, the crux of the problem is really the faverage that should be performed before the absolute value is taken.

[agramfort]

would you have time to open a PR so we can see a clear diff and look at a test?

…

[christian-oreilly]

Not in a near future but I can have a look at it when possible.

[Avoide]

Hi all,

I also had the same problem when working with resting-state epoched data.
The "error" described in this issue is not really an error but what @christian-oreilly describes is connectivity over time.
The current PLI and PLV implementation in spectral_connectivity() is connectivity over trials, which are suited for ERP data.

Mike X Cohen explains the difference quite nicely in his video here:
Youtube Link

I tried implementing my own version to calculate PLV and wPLI over time and you can see a simple working example here:

# Libraries
import numpy as np
import mne
from scipy import special
from mne.connectivity import spectral_connectivity

# Generate data
np.random.seed(42)
n_epochs = 5
n_channels = 3
n_times = 1000
data = np.random.rand(n_epochs, n_channels, n_times)

# Set sampling freq
sfreq = 250 # A reasonable random choice

# Choose what kind of data should be used for evaluation
# 0 = same trial repeated in all epochs
# 1 = 10Hz sinus waves with random phase differences in each channel and epoch
data_option = 1

if data_option == 0:
    # 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:
    # 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")

# Define freq bands
Freq_Bands = {"delta": [1.25, 4.0],
              "theta": [4.0, 8.0],
              "alpha": [8.0, 13.0],
              "beta": [13.0, 30.0],
              "gamma": [30.0, 49.0]}
n_freq_bands = len(Freq_Bands)
# Convert to tuples for the mne function
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))])

# Connectivity methods
connectivity_methods = ["plv","wpli"]
n_con_methods=len(connectivity_methods)

# Number of pairwise ch connections
n_ch_connections = scipy.special.comb(n_channels,2, exact=True, repetition=False)

# Pre-allocatate memory
con_data = np.zeros((n_con_methods,n_channels,n_channels,n_freq_bands))
con_data[con_data==0] = np.nan # nan matrix as 0 is meaningful

# Calculate PLV and wPLI - the MNE python implementation is over trials
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    data, method=connectivity_methods,
    mode="multitaper", sfreq=sfreq, fmin=fmin, fmax=fmax,
    faverage=True, verbose=0)
# Save the results in array
con_data[0,:,:,:] = con[0] # PLV
con_data[1,:,:,:] = con[1] # WPLI

print("Alpha PLV over trials")
print(con_data[0,:,:,2]) # all 1, since it is the same trial
# so PLV across trials should be 1

# PLV and wPLI across time
# Make linspace array for morlet waves
freq_centers = np.arange(fmin[0],fmax[-1]+0.25,0.25)
# Prepare Morlets
morlets = mne.time_frequency.tfr.morlet(sfreq,freq_centers,n_cycles=3)

# Make freqs array for indexing
freqs0 = [0]*n_freq_bands
for f in range(n_freq_bands):
    freqs0[f] = freq_centers[(freq_centers>=fmin[f]) & (freq_centers<=fmax[f])]

def calculate_PLV_WPLI_across_time(cwt_data):
    """

    Parameters
    ----------
    cwt_data : array, shape(n_channels, n_freq, n_times)
        The continuous wavelet transform of the data in one epoch

    Returns
    -------
    con_array : array, shape(2, n_channels, n_channels, n_freq_bands)
        The connectivity matrix. Only the lower diagonal is calculated.
        First axis indicates whether it is PLV (0) or wPLI(1)

    """
    n_ch, n_freq0, n_time0 = cwt_data.shape
    # The wavelet transform coefficients are complex numbers
    # The real part correspond to the amplitude and the imaginary part can be used to calculate the phase
    angles_data = np.apply_along_axis(np.angle,2,cwt_data)
    # Prepare array with phase differences between all combinations
    phase_diff_arr = np.zeros((n_ch_connections,n_freq0,n_time0))
    phase_diff_arr_ch_idx = [0]*n_ch_connections
    counter = 0
    for ch_c in range(n_ch):
        for ch_r in range(ch_c+1,n_ch): # only calculate lower diagonal
            phase_diff_arr[counter] = angles_data[ch_r]-angles_data[ch_c]
            phase_diff_arr_ch_idx[counter] = [ch_r,ch_c]
            counter += 1
            
    del angles_data # free up some memory
    # =========================================================================
    # PLV over time correspond to mean across time of the absolute value of
    # the circular length of the relative phases. So PLV will be 1 if
    # the phases of 2 signals maintain a constant lag
    # In equational form: PLV = 1/N * |sum(e^i(phase1-phase2))|
    # =========================================================================
    # Convert phase difference to complex part i(phase1-phase2)
    phase_diff_arr_im = 0*phase_diff_arr+1j*phase_diff_arr
    # Take the exponential
    expPhase_arr = np.apply_along_axis(np.exp,2,phase_diff_arr_im)
    # Take mean and then the absolute value
    meanexpPhase_arr = np.apply_along_axis(np.mean,2,expPhase_arr)
    PLV_arr = np.apply_along_axis(np.abs,1,meanexpPhase_arr)
    
    del phase_diff_arr_im, expPhase_arr # free up some memory
    # =========================================================================
    # PLI over time correspond to the sign of the sine of relative phase
    # differences. So PLI will be 1 if one signal is always leading or
    # lagging behind the other signal. But it is insensitive to changes in
    # relative phase, as long as it is the same signal that leads.
    # If 2 signals are almost in phase, they might shift between lead/lag
    # due to small fluctuations from noise. This would lead to unstable
    # estimation of "phase" synchronisation.
    # The wPLI tries to correct for this by weighting the PLI with the
    # magnitude of the lag, to attenuate noise sources giving rise to
    # near zero phase lag "synchronization"
    # In equational form: WPLI = |E{|phase_diff|*sign(phase_diff)}| / E{|phase_diff|}
    # =========================================================================
    # Calculate the magnitude of phase differences
    phase_diff_mag_arr = np.abs(np.sin(phase_diff_arr))
    # Calculate the signed phase difference (PLI)
    sign_phase_diff_arr = np.sign(np.sin(phase_diff_arr))
    # Calculate the nominator (abs and average across time)
    WPLI_nominator = np.abs(np.mean(phase_diff_mag_arr*sign_phase_diff_arr,axis=2))
    # Calculate denominator for normalization
    WPLI_denom = np.mean(phase_diff_mag_arr, axis=2)
    # Calculate WPLI
    WPLI_arr = WPLI_nominator/WPLI_denom
    
    del phase_diff_mag_arr, sign_phase_diff_arr, phase_diff_arr # free up some memory
    # Calculate mean for each freq band
    con_array0 = np.zeros((2,n_ch_connections,n_freq_bands))
    for f in range(n_freq_bands):
        freq_of_interest = freqs0[f] # find freqs in the freq band of interest
        freq_idx = [i in freq_of_interest for i in freq_centers] # get idx
        con_array0[0,:,f] = np.apply_along_axis(np.mean,1,PLV_arr[:,freq_idx])
        con_array0[1,:,f] = np.apply_along_axis(np.mean,1,WPLI_arr[:,freq_idx])
    
    # Save to final array with ch-ch connectivity in matrix form
    con_array = np.zeros((2,n_ch,n_ch,n_freq_bands))
    for com in range(n_ch_connections):
        ch_r = phase_diff_arr_ch_idx[com][0] # get row idx
        ch_c = phase_diff_arr_ch_idx[com][1] # get col idx
        con_array[0,ch_r,ch_c,:] = con_array0[0,com]
        con_array[1,ch_r,ch_c,:] = con_array0[1,com]

    return con_array

# Pre-allocate memory
con_data_time = np.zeros((n_con_methods,n_channels,n_channels,n_freq_bands))
con_data_time[con_data_time==0] = np.nan # nan matrix as 0 is meaningful

con_data1 = np.zeros((n_con_methods,n_epochs,n_channels,n_channels,n_freq_bands))

for epoch in range(n_epochs):
    # First the data in each epoch is retrieved
    temp_data = data[i]
    # Then continuous wavelet transform is used to decompose in time frequencies
    temp_data_cwt = mne.time_frequency.tfr.cwt(temp_data,morlets)
    # PLV and WPLI value is calculated across timepoints in each freq band and averaged into the 5 defined freq bands
    PLV_WPLI_con = calculate_PLV_WPLI_across_time(temp_data_cwt)
    # Save results
    con_data1[0,epoch,:,:,:] = PLV_WPLI_con[0] # phase locking value
    con_data1[1,epoch,:,:,:] = PLV_WPLI_con[1] # weighted phase lag index

# Take average across epochs for PLV and wPLI
con_data2 = np.mean(con_data1,axis=1)
# Save to final array
con_data_time[0,:,:,:] = con_data2[0] # phase locking value
con_data_time[1,:,:,:] = con_data2[1] # weighted phase lag index

print("Alpha PLV over time")
print(con_data_time[0,:,:,2])

Notice I added the option to choose what kind of data should be used. If data_option = 0 then the expected output is that PLV across trials are 1 since I just repeated the first epoch in all other epochs, while if data_option = 1 the PLV across time should be 1 as I use 10Hz sine waves with different phase differences.

It is still quite rough and not optimized but should give an idea of the procedure.
I am not completely certain about my implementation but the wPLI by following equation 1 in [1] and weighing by the absolute magnitude of the phase differences as described by [2]. PLV was calculated according to the equation in [3]. According to mne-tools/mne-python#8305 this equation should be equivalent to what is being used in spectral_connectivity() already, but my implementation over time gives another PLV than spectral_connectivity().

It would be helpful if anyone could have a look.

[1] Hardmeier, Martin, Florian Hatz, Habib Bousleiman, Christian Schindler, Cornelis Jan Stam, and Peter Fuhr. 2014. “Reproducibility of Functional Connectivity and Graph Measures Based on the Phase Lag Index (PLI) and Weighted Phase Lag Index (WPLI) Derived from High Resolution EEG.” PLoS ONE 9 (10). https://doi.org/10.1371/journal.pone.0108648.
[2] Vinck, Martin, Robert Oostenveld, Marijn Van Wingerden, Franscesco Battaglia, and Cyriel M.A. Pennartz. 2011. “An Improved Index of Phase-Synchronization for Electrophysiological Data in the Presence of Volume-Conduction, Noise and Sample-Size Bias.” NeuroImage 55 (4): 1548–65. https://doi.org/10.1016/j.neuroimage.2011.01.055.
[3] Lachaux, Jean Philippe, Eugenio Rodriguez, Jacques Martinerie, and Francisco J. Varela. 1999. “Measuring Phase Synchrony in Brain Signals.” Human Brain Mapping 8 (4): 194–208. https://doi.org/10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C.

[Avoide]

would you have time to open a PR so we can see a clear diff and look at a test?
…

@agramfort The code I wrote includes testing on some simple simulated data and an implementation, but as I wrote it needs to be optimized better and there is nothing wrong with the current spectral_connectivity(), but it is just an option that could be added to address connectivity over time.

[agramfort]

feel free to contribute this to mne when you find the time

…

[juancamiload]

Perhaps the implementation in HyPyP may help?
e.g., https://github.com/ppsp-team/HyPyP/blob/a7f78ef9886ca0621b0d2b5d5650eb1d8c903fc8/hypyp/analyses.py#L398

[adam2392]

Hi @Avoide I think the latest commit on main might have this resolved with the PR #104 and #115. WDYT?

mne-connectivity/mne_connectivity/spectral/time.py

Lines 202 to 204 in ed40aef

	'pli' : Phase Lag Index (PLI) :footcite:`StamEtAl2007` given by::
	PLI = |E[sign(Im(Sxy))]|

[Avoide]

Hi @Avoide I think the latest commit on main might have this resolved with the PR #104 and #115. WDYT?

mne-connectivity/mne_connectivity/spectral/time.py

Lines 202 to 204 in ed40aef

	'pli' : Phase Lag Index (PLI) :footcite:`StamEtAl2007` given by::
	PLI = |E[sign(Im(Sxy))]|

Hi @adam2392. I had a look at the PRs and also tested the functions from the newest (dev) version of mne-connectivity on the simulated toy-examples from my earlier response in this Issue, and got the expected output! (Low connectivity over time, around 0.1-0.2 PLV/wPLI for random data, while more than 0.97 PLV/wPLI for a sinus wave), so I agree this have now been resolved.

Great work @ruuskas! It was also much easier to use the functions compared to when I last tested them, with the new API. The only minor discrepancy I noticed was that in spectral_connectivity_epochs, the freqs are automatically calculated when I input fmin and fmax, while in spectral_connectivity_time you are required to specify freqs, despite giving fmin and fmax.

[ruuskas]

The only minor discrepancy I noticed was that in spectral_connectivity_epochs, the freqs are automatically calculated when I input fmin and fmax, while in spectral_connectivity_time you are required to specify freqs, despite giving fmin and fmax.

Hi @Avoide! Good to hear that it's also working for you.
The reason spectral_connectivity_time requires freqs is that setting them automatically is quite arbitrary and could potentially even lead to misinterpretation of the results. I would say that a 1 Hz grid is a good starting point.

This is different to spectral_connectivity_epochs, where freqs are just the output frequencies of SciPy's rfft function.

[adam2392]

Yeah good point, we should probably consolidate the API to require freqs for both. WDYT @larsoner ? I think we could still support freqs = 'auto' but perhaps now freqs should be a required parameter?

[ruuskas]

This is just my take, but as far as I understand that wouldn't be possible or practical.

Freqs is already required with spectral_connectivity_epochs in the Morlet wavelet mode. In the multitaper and Fourier modes, the spectrum is computed using scipy.fft.rfft, which does not allow specifying freqs.

[adam2392]

I'm going to close this issue then as resolved via #115. If there's more API discussion to be had regarding spectral_connectivity_epochs, we can open another issue.