From 7d29201e40d830e7cca644e89763b23335544d4c Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Mon, 7 Feb 2022 12:34:10 -0500 Subject: [PATCH 01/16] Add and document DPLIEst_ class in epochs.py --- doc/references.bib | 13 +++++++++++ mne_connectivity/spectral/epochs.py | 36 +++++++++++++++++++++++++++-- 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/doc/references.bib b/doc/references.bib index d889b6ea..24b48327 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -226,3 +226,16 @@ @article{ColcloughEtAl2015 year = {2015}, pages = {439--448} } + + @article{StamEtAl2012, + title={Go with the flow: Use of a directed phase lag index (dPLI) to characterize patterns of phase relations in a large-scale model of brain dynamics}, + volume={62}, + ISSN={1053-8119}, + DOI={10.1016/j.neuroimage.2012.05.050}, + number={3}, + journal={NeuroImage}, + author={Stam, C. J. and van Straaten, E. C. W.}, + year={2012}, + month={Sep}, + pages={1415–1428} +} \ No newline at end of file diff --git a/mne_connectivity/spectral/epochs.py b/mne_connectivity/spectral/epochs.py index 558a7c0b..4ce9d001 100644 --- a/mne_connectivity/spectral/epochs.py +++ b/mne_connectivity/spectral/epochs.py @@ -374,6 +374,31 @@ def compute_con(self, con_idx, n_epochs): self.con_scores[con_idx] = con +class _DPLIEst(_EpochMeanConEstBase): + """DPLI Estimator.""" + + name = 'DPLI' + + def __init__(self, n_cons, n_freqs, n_times): + super(_DPLIEst, self).__init__(n_cons, n_freqs, n_times) + + # allocate accumulator + self._acc = np.zeros(self.csd_shape) + + def accumulate(self, con_idx, csd_xy): + """Accumulate some connections.""" + self._acc[con_idx] += np.heaviside(np.imag(csd_xy), 0.5) + + def compute_con(self, con_idx, n_epochs): + """Compute final con. score for some connections.""" + if self.con_scores is None: + self.con_scores = np.zeros(self.csd_shape) + + con = self._acc[con_idx] / n_epochs + + self.con_scores[con_idx] = con + + class _WPLIEst(_EpochMeanConEstBase): """WPLI Estimator.""" @@ -687,7 +712,8 @@ def _get_and_verify_data_sizes(data, sfreq, n_signals=None, n_times=None, _CON_METHOD_MAP = {'coh': _CohEst, 'cohy': _CohyEst, 'imcoh': _ImCohEst, 'plv': _PLVEst, 'ciplv': _ciPLVEst, 'ppc': _PPCEst, 'pli': _PLIEst, 'pli2_unbiased': _PLIUnbiasedEst, - 'wpli': _WPLIEst, 'wpli2_debiased': _WPLIDebiasedEst} + 'dpli': _DPLIEst, 'wpli': _WPLIEst, + 'wpli2_debiased': _WPLIDebiasedEst} def _check_estimators(method, mode): @@ -749,7 +775,8 @@ def spectral_connectivity_epochs(data, names=None, method='coh', indices=None, %(names)s method : str | list of str Connectivity measure(s) to compute. These can be ``['coh', 'cohy', - 'imcoh', 'plv', 'ciplv', 'ppc', 'pli', 'wpli', 'wpli2_debiased']``. + 'imcoh', 'plv', 'ciplv', 'ppc', 'pli', 'dpli', 'wpli', + 'wpli2_debiased']``. indices : tuple of array | None Two arrays with indices of connections for which to compute connectivity. If None, all connections are computed. @@ -904,6 +931,11 @@ def spectral_connectivity_epochs(data, names=None, method='coh', indices=None, 'pli2_unbiased' : Unbiased estimator of squared PLI :footcite:`VinckEtAl2011`. + 'dpli' : Directed Phase Lag Index (DPLI) :footcite`StamEtAl2012` + given by:: + + DPLI = E[H(Im(Sxy))] + 'wpli' : Weighted Phase Lag Index (WPLI) :footcite:`VinckEtAl2011` given by:: From 0a7a9bd7ab9e6ffbfb5babf2bc9c3daf7861bc65 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Mon, 7 Feb 2022 12:35:11 -0500 Subject: [PATCH 02/16] Add DPLI parametrization to test_spectral_connectivity --- mne_connectivity/spectral/tests/test_spectral.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne_connectivity/spectral/tests/test_spectral.py b/mne_connectivity/spectral/tests/test_spectral.py index 0cdb698c..4aaa927a 100644 --- a/mne_connectivity/spectral/tests/test_spectral.py +++ b/mne_connectivity/spectral/tests/test_spectral.py @@ -151,7 +151,7 @@ def test_spectral_connectivity_parallel(method, mode, tmp_path): @pytest.mark.parametrize('method', ['coh', 'cohy', 'imcoh', 'plv', ['ciplv', 'ppc', 'pli', 'pli2_unbiased', - 'wpli', 'wpli2_debiased', 'coh']]) + 'dpli', 'wpli', 'wpli2_debiased', 'coh']]) @pytest.mark.parametrize('mode', ['multitaper', 'fourier', 'cwt_morlet']) def test_spectral_connectivity(method, mode): """Test frequency-domain connectivity methods.""" From 33b5bb9f7525f645daa1e0ca650356a828b3f589 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Thu, 10 Feb 2022 15:40:28 -0500 Subject: [PATCH 03/16] Fix citation --- mne_connectivity/spectral/epochs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mne_connectivity/spectral/epochs.py b/mne_connectivity/spectral/epochs.py index 4ce9d001..102d8cda 100644 --- a/mne_connectivity/spectral/epochs.py +++ b/mne_connectivity/spectral/epochs.py @@ -931,8 +931,8 @@ def spectral_connectivity_epochs(data, names=None, method='coh', indices=None, 'pli2_unbiased' : Unbiased estimator of squared PLI :footcite:`VinckEtAl2011`. - 'dpli' : Directed Phase Lag Index (DPLI) :footcite`StamEtAl2012` - given by:: + 'dpli' : Directed Phase Lag Index (DPLI) :footcite:`StamEtAl2012` + given by (where H is the Heaviside function):: DPLI = E[H(Im(Sxy))] From 0636caace7874152b94df7715b98122935bcba4a Mon Sep 17 00:00:00 2001 From: Adam Li Date: Thu, 10 Feb 2022 18:06:46 -0500 Subject: [PATCH 04/16] Fix circleCI --- mne_connectivity/spectral/time.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index 33dc4b37..1d9649c5 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -58,9 +58,9 @@ def spectral_connectivity_time(data, names=None, method='coh', indices=None, default, a 500ms smoothing is used. sm_freqs : int | 1 Number of points for frequency smoothing. By default, 1 is used which - is equivalent to no smoothing - kernel : {'square', 'hanning'} - Kernel type to use. Choose either 'square' or 'hanning' + is equivalent to no smoothing. + sm_kernel : {'square', 'hanning'} + Kernel type to use. Choose either 'square' or 'hanning'. mode : str, optional Spectrum estimation mode can be either: 'multitaper', or 'cwt_morlet'. From bf5a7032f90d59d7ed7efb3cdba92bb505260e7b Mon Sep 17 00:00:00 2001 From: Adam Li Date: Thu, 10 Feb 2022 21:47:02 -0500 Subject: [PATCH 05/16] Adding fix --- mne_connectivity/spectral/time.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index 1d9649c5..bdda44bb 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -58,9 +58,9 @@ def spectral_connectivity_time(data, names=None, method='coh', indices=None, default, a 500ms smoothing is used. sm_freqs : int | 1 Number of points for frequency smoothing. By default, 1 is used which - is equivalent to no smoothing. - sm_kernel : {'square', 'hanning'} - Kernel type to use. Choose either 'square' or 'hanning'. + is equivalent to no smoothing + kernel : {'square', 'hanning'} + Kernel type to use. Choose either 'square' or 'hanning' mode : str, optional Spectrum estimation mode can be either: 'multitaper', or 'cwt_morlet'. @@ -215,8 +215,18 @@ def spectral_connectivity_time(data, names=None, method='coh', indices=None, # compute time-resolved spectral connectivity conn_tr = _spectral_connectivity(data[epoch_idx, ...], **call_params) + print(data.shape) + print(conn.shape) + print([x.shape for x in conn_tr]) + + # average across tapers if necessary + conn_tr = np.stack(conn_tr, axis=1) + if all(x.ndim == 4 for x in conn_tr): + conn_tr = np.mean(conn, axis=1) + + print(conn_tr.shape) # merge results - conn[epoch_idx, ...] = np.stack(conn_tr, axis=1) + conn[epoch_idx, ...] = conn_tr # create a Connectivity container indices = 'symmetric' @@ -255,7 +265,7 @@ def _spectral_connectivity(data, method, kernel, foi_idx, out += [tfr_array_multitaper( data, sfreq, [f_c], n_cycles=float(n_c), time_bandwidth=mt, output='complex', decim=decim, n_jobs=n_jobs, **kw_mt)] - out = np.stack(out, axis=2).squeeze() + out = np.stack(out, axis=3).squeeze() elif isinstance(mt_bandwidth, (type(None), int, float)): out = tfr_array_multitaper( data, sfreq, freqs, n_cycles=n_cycles, @@ -264,6 +274,10 @@ def _spectral_connectivity(data, method, kernel, foi_idx, # get the supported connectivity function conn_func = {'coh': _coh, 'plv': _plv, 'sxy': _cs}[method] + + # XXX: averages over tapers if necessary, but I don't think this is correct + if out.ndim == 5: + out = np.mean(out, axis=2) # computes conn across trials this_conn = conn_func(out, kernel, foi_idx, source_idx, target_idx, From 77f8cc6241efe2de791b0061e51159f8f0c65ee0 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Thu, 10 Feb 2022 21:48:06 -0500 Subject: [PATCH 06/16] Fixed pep --- mne_connectivity/spectral/time.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index bdda44bb..988e32da 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -274,7 +274,7 @@ def _spectral_connectivity(data, method, kernel, foi_idx, # get the supported connectivity function conn_func = {'coh': _coh, 'plv': _plv, 'sxy': _cs}[method] - + # XXX: averages over tapers if necessary, but I don't think this is correct if out.ndim == 5: out = np.mean(out, axis=2) From 9c2f68882f24fbf2ad1f1995e21ff21c1672e057 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Fri, 11 Feb 2022 10:50:16 -0500 Subject: [PATCH 07/16] Fix doc for sm_kernel and sm_freqs in spectral_connectivity_time --- mne_connectivity/spectral/time.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index 988e32da..0ecf17ed 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -58,8 +58,8 @@ def spectral_connectivity_time(data, names=None, method='coh', indices=None, default, a 500ms smoothing is used. sm_freqs : int | 1 Number of points for frequency smoothing. By default, 1 is used which - is equivalent to no smoothing - kernel : {'square', 'hanning'} + is equivalent to no smoothing. + sm_kernel : {'square', 'hanning'} Kernel type to use. Choose either 'square' or 'hanning' mode : str, optional Spectrum estimation mode can be either: 'multitaper', or From c3ba72f8e2cfe283f5ae61a81ac79f100c85b566 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Thu, 3 Mar 2022 21:49:55 -0500 Subject: [PATCH 08/16] Add wPLI/dPLI/PLI example --- examples/dpli_wpli_pli.py | 381 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 381 insertions(+) create mode 100644 examples/dpli_wpli_pli.py diff --git a/examples/dpli_wpli_pli.py b/examples/dpli_wpli_pli.py new file mode 100644 index 00000000..8b88ab63 --- /dev/null +++ b/examples/dpli_wpli_pli.py @@ -0,0 +1,381 @@ +''' +============================= +Comparing PLI, wPLI, and dPLI +============================= + +This example demonstrates the different connectivity information captured by +the phase lag index (PLI) :footcite:`StamEtAl2007`, weighted phase lag index +(wPLI) :footcite:`VinckEtAl2011`, and directed phase lag index (dPLI) +:footcite:`StamEtAl2012` on simulated data. +''' + +# Authors: Kenji Marshall +# Charlotte Maschke +# Stefanie Blain-Moraes +# +# License: BSD (3-clause) + + +import mne +import numpy as np +import matplotlib.pyplot as plt + +from mne_connectivity import spectral_connectivity_epochs +from mne.datasets import sample + +############################################################################### +# Background +# ---------- +# +# The formulae for PLI, wPLI, and dPLI are given below. In these equations, +# :math:`X_{ij}` is the cross-spectral density (CSD) between two signals +# :math:`i, j`. Importantly, the imaginary component of the CSD is maximal +# when the signals have a phase difference given by :math:`k\pi+\frac{\pi}{2}`, +# and is :math:`0` when the phase difference is given by :math:`k\pi` (where +# :math:`k \in \mathbb{Z}`). This property provides protection against +# recognizing volume conduction effects as connectivity, and is the backbone +# for these methods :footcite:`VinckEtAl2011`. +# +# :math:`PLI = |E[sgn(\mathcal{I}(X_{ij}))]|` :footcite:`StamEtAl2007` +# +# :math:`wPLI = \frac{|E[\mathcal{I}(X_{ij})]|}{E[|\mathcal{I}(X_{ij})|]}` +# :footcite:`VinckEtAl2011` +# +# :math:`dPLI = E[\mathcal{H}(\mathcal{I}(X_{ij}))]` :footcite:`StamEtAl2012` +# +# All three of these metrics are bounded in the range :math:`[0, 1]`. +# +# * For PLI, :math:`0` means that signal :math:`i` leads and lags signal +# :math:`j` equally often, while a value greater than :math:`0` means that +# there is an imbalance in the likelihood for signal :math:`i` to be leading +# or lagging. A value of :math:`1` means that signal :math:`i` only leads or +# only lags signal :math:`j`. +# * For wPLI, :math:`0` means that the total weight (not the quantity) of all +# leading relationships equals the total weight of lagging relationships, +# while a value greater than :math:`0` means that there is an imbalance +# between these weights. A value of :math:`1`, just as in PLI, means that +# signal :math:`i` only leads or only lags signal :math:`j`. +# * With dPLI, we gain the ability to distinguish whether signal :math:`i` is +# leading or lagging signal :math:`j`, complementing the information provided +# by PLI or wPLI. A value of :math:`0.5` means that signal :math:`i` leads +# and lags signal :math:`j` equally often. A value in the range +# :math:`(0.5, 1.0]` means that signal :math:`i` leads signal :math:`j` more +# often than it lags, with a value of :math:`1` meaning that signal :math:`i` +# always leads signal :math:`j`. A value in the range :math:`[0.0, 0.5)` +# means that signal :math:`i` lags signal :math:`j` more often than it leads, +# with a value of :math:`0` meaning that signal :math:`i` always lags signal +# :math:`j`. The PLI can actually be extracted from the dPLI by the +# relationship :math:`PLI = 2|dPLI - 0.5|`, but this relationship is not +# invertible (dPLI can not be estimated from the PLI). +# +# +# Overall, these different approaches are closely related but have subtle +# differences, as will be demonstrated throughout the rest of this example. + +############################################################################### +# Capturing Leading/Lagging Phase Relationships with dPLI +# ------------------------------------------------------- +# +# The main advantage of dPLI is that it's *directed*, meaning it can +# differentiate between phase relationships which are leading or lagging. +# To illustrate this, we generate sinusoids with Gaussian noise. In particular, +# we generate signals with phase differences of +# :math:`[-\pi, -\frac{\pi}{2}, 0, \frac{\pi}{2}, \pi]` relative to a reference +# signal. A negative difference means that the reference signal is lagging the +# other signal. + + +fs = 250 # sampling rate (Hz) +n_e = 300 # number of epochs +T = 10 # length of epochs (s) +f = 10 # frequency of sinusoids (Hz) +t = np.arange(0, T, 1 / fs) +A = 1 # noise amplitude +sigma = 0.5 # Gaussian noise variance + +data = [] + +phase_differences = [0, -np.pi, -np.pi / 2, 0, np.pi / 2, np.pi] +for ps in zip(phase_differences): + sig = [] + for _ in range(n_e): + sig.append(np.sin(2 * np.pi * f * t - ps) + + A * np.random.normal(0, sigma, size=t.shape)) + data.append(sig) + +data = np.swapaxes(np.array(data), 0, 1) # make epochs the first dimension + +############################################################################### +# A snippet of this simulated data is shown below. The blue signal is the +# reference signal. + +# %% +fig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True) +ax[0].plot(t[:fs], data[0, 0, :fs], label="Reference") +ax[0].plot(t[:fs], data[0, 2, :fs]) + +ax[0].set_title(r"Phase Lagging ($-\pi/2$ Phase Difference)") +ax[0].set_xlabel("Time (s)") +ax[0].set_ylabel("Signal") +ax[0].legend() + +ax[1].plot(t[:fs], data[0, 0, :fs], label="Reference") +ax[1].plot(t[:fs], data[0, 4, :fs]) +ax[1].set_title(r"Phase Leading ($\pi/2$ Phase Difference)") +ax[1].set_xlabel("Time (s)") + +plt.show() + +############################################################################### +# We will now compute PLI, wPLI, and dPLI for each phase relationship. + +# %% +conn = [] +indices = ([0, 0, 0, 0, 0], [1, 2, 3, 4, 5]) +for method in ['pli', 'wpli', 'dpli']: + conn.append( + spectral_connectivity_epochs( + data, method=method, sfreq=fs, indices=indices, + fmin=9, fmax=11, faverage=True).get_data()[:, 0]) +conn = np.array(conn) + +############################################################################### +# The estimated connectivites are shown in the figure below, which provides +# insight into the differences between PLI/wPLI, and dPLI. +# +# +# **Similarities Of All Measures** +# +# * Capture presence of connectivity in same situations (phase difference of +# :math:`\pm\frac{\pi}{2}`) +# * Do not predict connectivity when phase difference is a multiple of +# :math:`\pi` +# * Bounded between :math:`0` and :math:`1` +# +# **How dPLI is Different Than PLI/wPLI** +# +# * Null connectivity is :math:`0` for PLI and wPLI, but :math:`0.5` for dPLI +# * dPLI differentiates whether the reference signal is leading or lagging the +# other signal (lagging if :math:`0 <= dPlI < 0.5`, leading if +# :math:`0.5 < dPLI <= 1.0`) + + +x = np.arange(5) + +plt.figure() +plt.bar(x - 0.2, conn[0], 0.2, align='center', label="PLI") +plt.bar(x, conn[1], 0.2, align='center', label="wPLI") +plt.bar(x + 0.2, conn[2], 0.2, align='center', label="dPLI") + +plt.title("Connectivity Estimation Comparison") +plt.xticks(x, (r"$-\pi$", r"$-\pi/2$", r"$0$", r"$\pi/2$", r"$\pi$")) +plt.legend() +plt.xlabel("Phase Difference") +plt.ylabel("Estimated Connectivity") + +plt.show() + +############################################################################### +# Robustness to Outliers and Noise with wPLI +# ------------------------------------------ +# +# The previous experiment illustrated the advantages conferred by dPLI when +# differentiating leading and lagging phase relationships. This experiment +# will now focus on understanding the advantages of wPLI, and explore how it +# extends upon PLI. +# +# The main difference between PLI and wPLI is in how different phase +# relationships are *weighted*. In PLI, phase differences are weighted as +# :math:`-1` or :math:`1` according to their sign. In wPLI, phase differences +# are weighted based on their value, meaning that phase differences closer to +# :math:`\pm\frac{\pi}{2}` are weighted more heavily than those close to +# :math:`0` or any other multiple of :math:`\pi`. +# +# This avoids a discontinuity at the transition between positive and negative +# phase, treating all phase differences near this transition in a similar way. +# This provides some robustness against outliers and noise when estimating +# connectivity. For instance, volume conduction can distort EEG/MEG recordings, +# wherein signals emanating from the same neural source will be picked up by +# multiple sensors on the scalp. This can effect connectivity estimations, +# bringing the relative phase differences between two signals close to +# :math:`0`. wPLI minimizes the contribution of phase relationships that are +# small but non-zero (and may thus be attributed to volume conduciton), while +# PLI weighs these in the same way as phase relationships of +# :math:`\pm\frac{\pi}{2}`. +# +# To demonstrate this, we recreate a result from (Vinck et al, 2011) +# :footcite:`VinckEtAl2011`. Two sinusoids are simulated, where the phase +# difference for half of the epochs is :math:`\frac{\pi}{2}`, and is +# :math:`-\frac{\pi}{100}` for the others. We also explore the effect of +# applying uniform noise to this phase difference. + +# %% +n_noise = 41 # amount of noise amplitude samples in [0, 4] +data = [[]] + +# Generate reference +for _ in range(n_e): + data[0].append(np.sin(2 * np.pi * f * t)) + +A_list = np.linspace(0, 4, n_noise) + +for A in A_list: + sig = [] + # Generate other signal + for _ in range(int(n_e / 2)): # phase difference -pi/100 + sig.append(np.sin(2 * np.pi * f * t + np.pi / + 100 + A * np.random.uniform(-1, 1))) + for _ in range(int(n_e / 2), n_e): # phase difference pi/2 + sig.append(np.sin(2 * np.pi * f * t - np.pi / + 2 + A * np.random.uniform(-1, 1))) + data.append(sig) + +data = np.swapaxes(np.array(data), 0, 1) + +# Visualize the data +fig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True) +ax[0].plot(t[:10], data[0, 0, :10], label="Reference") +ax[0].plot(t[:10], data[0, 1, :10]) + +ax[0].set_title(r"Phase Lagging ($-\pi/100$ Phase Difference)") +ax[0].set_xlabel("Time (s)") +ax[0].set_ylabel("Signal") +ax[0].legend() + +ax[1].plot(t[:fs], data[0, 0, :fs], label="Reference") +ax[1].plot(t[:fs], data[0, -1, :fs]) +ax[1].set_title(r"Phase Leading ($\pi/2$ Phase Difference)") +ax[1].set_xlabel("Time (s)") + +plt.show() + +############################################################################### +# We can now compute PLI and wPLI + +conn = [] +indices = ([0] * n_noise, np.arange(1, n_noise + 1)) +for method in ['pli', 'wpli']: + conn.append( + spectral_connectivity_epochs( + data, method=method, sfreq=fs, indices=indices, + fmin=9, fmax=11, faverage=True).get_data()[:, 0]) +conn = np.array(conn) + +############################################################################### +# The results from the simulation are shown in the figure below. In the case +# without noise, the difference between wPLI and PLI is made obvious. In PLI, +# no connectivity is detected, as the :math:`-\frac{\pi}{100}` phase +# differences are weighted in the exact same way as the :math:`\frac{\pi}{2}` +# relationships. wPLI is able to avoid the cancellation of the +# :math:`\frac{\pi}{2}` relationships. +# +# As noise gets added, PLI increases since the :math:`-\frac{\pi}{100}` +# relationships are made positive more often than the :math:`\frac{\pi}{2}` +# relationships are made negative. However, wPLI maintains an advantage in +# its ability to distinguish the underlying structure. Beyond a certain point, +# the noise dominates any pre-defined structure, and both methods behave +# similarly, tending toward :math:`0`. For a more detailed analysis of this +# result and the properties of wPLI, please refer to (Vinck et al, 2011) +# :footcite:`VinckEtAl2011`. + +plt.figure() +plt.plot(A_list, conn[0], "o-", label="PLI") +plt.plot(A_list, conn[1], "o-", label="wPLI") +plt.legend() +plt.xlabel("Noise Amplitude") +plt.ylabel("Connectivity Measure") +plt.title("wPLI and PLI Under Increasing Noise") +plt.show() + +############################################################################### +# Demo On MEG Data +# ---------------- +# +# To finish this example, we also quickly demonstrate these methods on some +# sample MEG data recorded during visual stimulation. + +data_path = sample.data_path() +raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' +event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' +raw = mne.io.read_raw_fif(raw_fname) +events = mne.read_events(event_fname) + + +# Select gradiometers +picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True, + exclude='bads') + +# Create epochs +event_id, tmin, tmax = 3, -0.2, 1.5 # need a long enough epoch for 5 cycles +epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, + baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6)) +epochs.load_data().pick_types(meg='grad') # just keep MEG and no EOG now + +fmin, fmax = 4., 9. # compute connectivity within 4-9 Hz +sfreq = raw.info['sfreq'] # the sampling frequency +tmin = 0.0 # exclude the baseline period + +# Compute PLI, wPLI, and dPLI +con_pli = spectral_connectivity_epochs( + epochs, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin, + fmax=fmax, faverage=True, tmin=tmin, mt_adaptive=False, n_jobs=1) + +con_wpli = spectral_connectivity_epochs( + epochs, method='wpli', mode='multitaper', sfreq=sfreq, fmin=fmin, + fmax=fmax, faverage=True, tmin=tmin, mt_adaptive=False, n_jobs=1) + +con_dpli = spectral_connectivity_epochs( + epochs, method='dpli', mode='multitaper', sfreq=sfreq, fmin=fmin, + fmax=fmax, faverage=True, tmin=tmin, mt_adaptive=False, n_jobs=1) + +############################################################################### +# In this example, there is strong connectivity between sensors 190-200 and +# sensors 110-160. +# +# Moreover, after observing the presence of connectivity, dPLI can be used to +# ascertain the direction of the phase relationship. Here, it appears that the +# dPLI connectivity in this area is less than :math:`0.5`, and thus sensors +# 190-200 are lagging sensors 110-160. +# +# In keeping with the previous simulation, we can see that wPLI identifies +# stronger connectivity relationships than PLI. This is due to its robustness +# against volume conduction effects decreasing the detected connectivity +# strength, as was mentioned earlier. + +fig, axs = plt.subplots(1, 3, figsize=(14, 5), sharey=True) +axs[0].imshow(con_pli.get_data('dense'), vmin=0, vmax=1) +axs[0].set_title("PLI") +axs[0].set_ylabel("Sensor 1") +axs[0].set_xlabel("Sensor 2") + +axs[1].imshow(con_wpli.get_data('dense'), vmin=0, vmax=1) +axs[1].set_title("wPLI") +axs[1].set_xlabel("Sensor 2") + +im = axs[2].imshow(con_dpli.get_data('dense'), vmin=0, vmax=1) +axs[2].set_title("dPLI") +axs[2].set_xlabel("Sensor 2") + +fig.colorbar(im, ax=axs.ravel()) +plt.show() + +############################################################################### +# Conclusions +# ----------- +# +# Both wPLI and dPLI are extensions upon the original PLI method, and provide +# complementary information about underlying connectivity. +# +# * To identify the presence of an underlying phase relationship, wPLI is the +# method of choice for most researchers as it provides an improvement in +# robustness over the original PLI method +# * To know the directionality of the connectivity identified by wPLI, dPLI +# should be used +# +# Ultimately, these methods work great together, providing a comprehensive +# estimate of phase-based connectivity. + +############################################################################### +# References +# ---------- +# .. footbibliography:: From ae12528debefcdbf0bbdf8c964a87f9c8a1a58fe Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Thu, 3 Mar 2022 21:50:04 -0500 Subject: [PATCH 09/16] Update whatsnew --- doc/authors.inc | 3 ++- doc/whats_new.rst | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/doc/authors.inc b/doc/authors.inc index 255360e2..a92c90e0 100644 --- a/doc/authors.inc +++ b/doc/authors.inc @@ -1,4 +1,5 @@ .. _Adam Li: https://github.com/adam2392 .. _Eric Larson: https://github.com/larsoner .. _Britta Westner: https://github.com/britta-wstnr -.. _Alexander Kroner: https://github.com/alexanderkroner \ No newline at end of file +.. _Alexander Kroner: https://github.com/alexanderkroner +.. _Kenji Marshall: https://github.com/kenjimarshall \ No newline at end of file diff --git a/doc/whats_new.rst b/doc/whats_new.rst index dd3a277c..b05c4c07 100644 --- a/doc/whats_new.rst +++ b/doc/whats_new.rst @@ -29,6 +29,7 @@ Enhancements - Added :func:`mne_connectivity.select_order` for helping to select VAR order using information criterion, by `Adam Li`_ (:gh:`46`) - All connectivity functions retain ``events``, ``event_id`` and ``metadata`` from `mne.Epochs` objects as input and is stored as part of the connectivity object, by `Adam Li`_ (:gh:`58`) - Add spectral connectivity over time function :func:`mne_connectivity.spectral_connectivity_time`, by `Adam Li`_ (:gh:`67`) +- Add directed phase lag index (dPLI) as a method in :func:`mne_connectivity.spectral_connectivity_epochs` with a corresponding example by `Kenji Marshall`_ (:gh:`79`) Bug ~~~ @@ -51,6 +52,7 @@ Authors * `Adam Li`_ * `Eric Larson`_ * `Alexander Kroner`_ +* `Kenji Marshall`_ :doc:`Find out what was new in previous releases ` From 64fdbf151bcc24226353ee25c93d207dd1393336 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Tue, 29 Mar 2022 13:29:24 -0400 Subject: [PATCH 10/16] Update mne_connectivity/spectral/time.py Co-authored-by: Adam Li --- mne_connectivity/spectral/time.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index d9b2b709..2e564abc 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -215,16 +215,6 @@ def spectral_connectivity_time(data, names=None, method='coh', indices=None, # compute time-resolved spectral connectivity conn_tr = _spectral_connectivity(data[epoch_idx, ...], **call_params) - print(data.shape) - print(conn.shape) - print([x.shape for x in conn_tr]) - - # average across tapers if necessary - conn_tr = np.stack(conn_tr, axis=1) - if all(x.ndim == 4 for x in conn_tr): - conn_tr = np.mean(conn, axis=1) - - print(conn_tr.shape) # merge results conn[epoch_idx, ...] = conn_tr From 8540aeee302399d87cec4641e4974c38222c4aa2 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Tue, 29 Mar 2022 13:29:32 -0400 Subject: [PATCH 11/16] Update mne_connectivity/spectral/time.py Co-authored-by: Adam Li --- mne_connectivity/spectral/time.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index 2e564abc..623de03e 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -216,7 +216,7 @@ def spectral_connectivity_time(data, names=None, method='coh', indices=None, conn_tr = _spectral_connectivity(data[epoch_idx, ...], **call_params) # merge results - conn[epoch_idx, ...] = conn_tr + conn[epoch_idx, ...] = np.stack(conn_tr, axis=1) # create a Connectivity container indices = 'symmetric' From 05e6fbb0590bfc8360805e0bf13a468af7942a65 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Tue, 29 Mar 2022 13:30:46 -0400 Subject: [PATCH 12/16] Add imaginary and heaviside definitions --- examples/dpli_wpli_pli.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/dpli_wpli_pli.py b/examples/dpli_wpli_pli.py index 8b88ab63..c30c4a7e 100644 --- a/examples/dpli_wpli_pli.py +++ b/examples/dpli_wpli_pli.py @@ -34,7 +34,9 @@ # and is :math:`0` when the phase difference is given by :math:`k\pi` (where # :math:`k \in \mathbb{Z}`). This property provides protection against # recognizing volume conduction effects as connectivity, and is the backbone -# for these methods :footcite:`VinckEtAl2011`. +# for these methods :footcite:`VinckEtAl2011`. In the equations below, +# :math:`\mathcal{I}` refers to the imaginary component, while +# :math:`\mathcal{H}` refers to the Heaviside step function. # # :math:`PLI = |E[sgn(\mathcal{I}(X_{ij}))]|` :footcite:`StamEtAl2007` # From 3ab927ba5457bfc0705a5131bd82b037623fbea8 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Tue, 29 Mar 2022 13:34:47 -0400 Subject: [PATCH 13/16] Define sgn --- examples/dpli_wpli_pli.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/dpli_wpli_pli.py b/examples/dpli_wpli_pli.py index c30c4a7e..04c18bad 100644 --- a/examples/dpli_wpli_pli.py +++ b/examples/dpli_wpli_pli.py @@ -35,8 +35,9 @@ # :math:`k \in \mathbb{Z}`). This property provides protection against # recognizing volume conduction effects as connectivity, and is the backbone # for these methods :footcite:`VinckEtAl2011`. In the equations below, -# :math:`\mathcal{I}` refers to the imaginary component, while -# :math:`\mathcal{H}` refers to the Heaviside step function. +# :math:`\mathcal{I}` refers to the imaginary component, +# :math:`\mathcal{H}` refers to the Heaviside step function, and +# :math:`sgn` refers to the sign function # # :math:`PLI = |E[sgn(\mathcal{I}(X_{ij}))]|` :footcite:`StamEtAl2007` # From eb7ecdb34ffb6c9d8239d57edf5a63f7bd7c5783 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 29 Mar 2022 14:00:05 -0400 Subject: [PATCH 14/16] Update mne_connectivity/spectral/time.py --- mne_connectivity/spectral/time.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index 623de03e..965779e7 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -269,10 +269,6 @@ def _spectral_connectivity(data, method, kernel, foi_idx, # get the supported connectivity function conn_func = {'coh': _coh, 'plv': _plv, 'sxy': _cs}[method] - # XXX: averages over tapers if necessary, but I don't think this is correct - if out.ndim == 5: - out = np.mean(out, axis=2) - # computes conn across trials # TODO: This is wrong -- it averages in the complex domain (over tapers). # What it *should* do is compute the conn for each taper, then average From 8c6ded72a2989a4f84272a9717f8865345352e26 Mon Sep 17 00:00:00 2001 From: Kenji Marshall Date: Tue, 29 Mar 2022 16:58:22 -0400 Subject: [PATCH 15/16] Fix lagging/leading plot --- examples/dpli_wpli_pli.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/dpli_wpli_pli.py b/examples/dpli_wpli_pli.py index 04c18bad..6865a346 100644 --- a/examples/dpli_wpli_pli.py +++ b/examples/dpli_wpli_pli.py @@ -37,7 +37,7 @@ # for these methods :footcite:`VinckEtAl2011`. In the equations below, # :math:`\mathcal{I}` refers to the imaginary component, # :math:`\mathcal{H}` refers to the Heaviside step function, and -# :math:`sgn` refers to the sign function +# :math:`sgn` refers to the sign function. # # :math:`PLI = |E[sgn(\mathcal{I}(X_{ij}))]|` :footcite:`StamEtAl2007` # @@ -238,7 +238,7 @@ # Visualize the data fig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True) ax[0].plot(t[:10], data[0, 0, :10], label="Reference") -ax[0].plot(t[:10], data[0, 1, :10]) +ax[0].plot(t[:10], data[1, 1, :10]) ax[0].set_title(r"Phase Lagging ($-\pi/100$ Phase Difference)") ax[0].set_xlabel("Time (s)") @@ -246,7 +246,7 @@ ax[0].legend() ax[1].plot(t[:fs], data[0, 0, :fs], label="Reference") -ax[1].plot(t[:fs], data[0, -1, :fs]) +ax[1].plot(t[:fs], data[-1, 1, :fs]) ax[1].set_title(r"Phase Leading ($\pi/2$ Phase Difference)") ax[1].set_xlabel("Time (s)") From 3d91c6eafdd967a212d5af524223f7346bf2229c Mon Sep 17 00:00:00 2001 From: Adam Li Date: Wed, 30 Mar 2022 11:22:48 -0400 Subject: [PATCH 16/16] Update doc/whats_new.rst --- doc/whats_new.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/whats_new.rst b/doc/whats_new.rst index b73251e0..8ec53449 100644 --- a/doc/whats_new.rst +++ b/doc/whats_new.rst @@ -41,6 +41,9 @@ Authors ~~~~~~~ * `Kenji Marshall`_ +* `Adam Li`_ +* `Alex Rockhill`_ +* `Szonja Weigl`_ :doc:`Find out what was new in previous releases `