{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Comparing PLI, wPLI, and dPLI\n\nThis example demonstrates the different connectivity information captured by\nthe phase lag index (PLI) :footcite:`StamEtAl2007`, weighted phase lag index\n(wPLI) :footcite:`VinckEtAl2011`, and directed phase lag index (dPLI)\n:footcite:`StamEtAl2012` on simulated data.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Authors: Kenji Marshall \n# Charlotte Maschke \n# Stefanie Blain-Moraes \n#\n# License: BSD (3-clause)\n\n\nimport matplotlib.pyplot as plt\nimport mne\nimport numpy as np\nfrom mne.datasets import sample\n\nfrom mne_connectivity import spectral_connectivity_epochs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Background\n\nThe formulae for PLI, wPLI, and dPLI are given below. In these equations,\n$X_{ij}$ is the cross-spectral density (CSD) between two signals\n$i, j$. Importantly, the imaginary component of the CSD is maximal\nwhen the signals have a phase difference given by $k\\pi+\\frac{\\pi}{2}$,\nand is $0$ when the phase difference is given by $k\\pi$ (where\n$k \\in \\mathbb{Z}$). This property provides protection against\nrecognizing volume conduction effects as connectivity, and is the backbone\nfor these methods :footcite:`VinckEtAl2011`. In the equations below,\n$\\mathcal{I}$ refers to the imaginary component,\n$\\mathcal{H}$ refers to the Heaviside step function, and\n$sgn$ refers to the sign function.\n\n$PLI = |E[sgn(\\mathcal{I}(X_{ij}))]|$ :footcite:`StamEtAl2007`\n\n$wPLI = \\frac{|E[\\mathcal{I}(X_{ij})]|}{E[|\\mathcal{I}(X_{ij})|]}$\n:footcite:`VinckEtAl2011`\n\n$dPLI = E[\\mathcal{H}(\\mathcal{I}(X_{ij}))]$ :footcite:`StamEtAl2012`\n\nAll three of these metrics are bounded in the range $[0, 1]$.\n\n* For PLI, $0$ means that signal $i$ leads and lags signal\n $j$ equally often, while a value greater than $0$ means that\n there is an imbalance in the likelihood for signal $i$ to be leading\n or lagging. A value of $1$ means that signal $i$ only leads or\n only lags signal $j$.\n* For wPLI, $0$ means that the total weight (not the quantity) of all\n leading relationships equals the total weight of lagging relationships,\n while a value greater than $0$ means that there is an imbalance\n between these weights. A value of $1$, just as in PLI, means that\n signal $i$ only leads or only lags signal $j$.\n* With dPLI, we gain the ability to distinguish whether signal $i$ is\n leading or lagging signal $j$, complementing the information provided\n by PLI or wPLI. A value of $0.5$ means that signal $i$ leads\n and lags signal $j$ equally often. A value in the range\n $(0.5, 1.0]$ means that signal $i$ leads signal $j$ more\n often than it lags, with a value of $1$ meaning that signal $i$\n always leads signal $j$. A value in the range $[0.0, 0.5)$\n means that signal $i$ lags signal $j$ more often than it leads,\n with a value of $0$ meaning that signal $i$ always lags signal\n $j$. The PLI can actually be extracted from the dPLI by the\n relationship $PLI = 2|dPLI - 0.5|$, but this relationship is not\n invertible (dPLI can not be estimated from the PLI).\n\n\nOverall, these different approaches are closely related but have subtle\ndifferences, as will be demonstrated throughout the rest of this example.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Capturing Leading/Lagging Phase Relationships with dPLI\n\nThe main advantage of dPLI is that it's *directed*, meaning it can\ndifferentiate between phase relationships which are leading or lagging.\nTo illustrate this, we generate sinusoids with Gaussian noise. In particular,\nwe generate signals with phase differences of\n$[-\\pi, -\\frac{\\pi}{2}, 0, \\frac{\\pi}{2}, \\pi]$ relative to a reference\nsignal. A negative difference means that the reference signal is lagging the\nother signal.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fs = 250 # sampling rate (Hz)\nn_e = 300 # number of epochs\nT = 10 # length of epochs (s)\nf = 10 # frequency of sinusoids (Hz)\nt = np.arange(0, T, 1 / fs)\nA = 1 # noise amplitude\nsigma = 0.5 # Gaussian noise variance\n\ndata = []\n\nphase_differences = [0, -np.pi, -np.pi / 2, 0, np.pi / 2, np.pi]\nfor ps in zip(phase_differences):\n sig = []\n for _ in range(n_e):\n sig.append(\n np.sin(2 * np.pi * f * t - ps)\n + A * np.random.normal(0, sigma, size=t.shape)\n )\n data.append(sig)\n\ndata = np.swapaxes(np.array(data), 0, 1) # make epochs the first dimension"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A snippet of this simulated data is shown below. The blue signal is the\nreference signal.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True)\nax[0].plot(t[:fs], data[0, 0, :fs], label=\"Reference\")\nax[0].plot(t[:fs], data[0, 2, :fs])\n\nax[0].set_title(r\"Phase Lagging ($-\\pi/2$ Phase Difference)\")\nax[0].set_xlabel(\"Time (s)\")\nax[0].set_ylabel(\"Signal\")\nax[0].legend()\n\nax[1].plot(t[:fs], data[0, 0, :fs], label=\"Reference\")\nax[1].plot(t[:fs], data[0, 4, :fs])\nax[1].set_title(r\"Phase Leading ($\\pi/2$ Phase Difference)\")\nax[1].set_xlabel(\"Time (s)\")\n\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will now compute PLI, wPLI, and dPLI for each phase relationship.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"conn = []\nindices = ([0, 0, 0, 0, 0], [1, 2, 3, 4, 5])\nfor method in [\"pli\", \"wpli\", \"dpli\"]:\n conn.append(\n spectral_connectivity_epochs(\n data,\n method=method,\n sfreq=fs,\n indices=indices,\n fmin=9,\n fmax=11,\n faverage=True,\n ).get_data()[:, 0]\n )\nconn = np.array(conn)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The estimated connectivites are shown in the figure below, which provides\ninsight into the differences between PLI/wPLI, and dPLI.\n\n\n**Similarities Of All Measures**\n\n* Capture presence of connectivity in same situations (phase difference of\n $\\pm\\frac{\\pi}{2}$)\n* Do not predict connectivity when phase difference is a multiple of\n $\\pi$\n* Bounded between $0$ and $1$\n\n**How dPLI is Different Than PLI/wPLI**\n\n* Null connectivity is $0$ for PLI and wPLI, but $0.5$ for dPLI\n* dPLI differentiates whether the reference signal is leading or lagging the\n other signal (lagging if $0 <= dPlI < 0.5$, leading if\n $0.5 < dPLI <= 1.0$)\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = np.arange(5)\n\nplt.figure()\nplt.bar(x - 0.2, conn[0], 0.2, align=\"center\", label=\"PLI\")\nplt.bar(x, conn[1], 0.2, align=\"center\", label=\"wPLI\")\nplt.bar(x + 0.2, conn[2], 0.2, align=\"center\", label=\"dPLI\")\n\nplt.title(\"Connectivity Estimation Comparison\")\nplt.xticks(x, (r\"$-\\pi$\", r\"$-\\pi/2$\", r\"$0$\", r\"$\\pi/2$\", r\"$\\pi$\"))\nplt.legend()\nplt.xlabel(\"Phase Difference\")\nplt.ylabel(\"Estimated Connectivity\")\n\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Robustness to Outliers and Noise with wPLI\n\nThe previous experiment illustrated the advantages conferred by dPLI when\ndifferentiating leading and lagging phase relationships. This experiment\nwill now focus on understanding the advantages of wPLI, and explore how it\nextends upon PLI.\n\nThe main difference between PLI and wPLI is in how different phase\nrelationships are *weighted*. In PLI, phase differences are weighted as\n$-1$ or $1$ according to their sign. In wPLI, phase differences\nare weighted based on their value, meaning that phase differences closer to\n$\\pm\\frac{\\pi}{2}$ are weighted more heavily than those close to\n$0$ or any other multiple of $\\pi$.\n\nThis avoids a discontinuity at the transition between positive and negative\nphase, treating all phase differences near this transition in a similar way.\nThis provides some robustness against outliers and noise when estimating\nconnectivity. For instance, volume conduction can distort EEG/MEG recordings,\nwherein signals emanating from the same neural source will be picked up by\nmultiple sensors on the scalp. This can effect connectivity estimations,\nbringing the relative phase differences between two signals close to\n$0$. wPLI minimizes the contribution of phase relationships that are\nsmall but non-zero (and may thus be attributed to volume conduciton), while\nPLI weighs these in the same way as phase relationships of\n$\\pm\\frac{\\pi}{2}$.\n\nTo demonstrate this, we recreate a result from (Vinck et al, 2011)\n:footcite:`VinckEtAl2011`. Two sinusoids are simulated, where the phase\ndifference for half of the epochs is $\\frac{\\pi}{2}$, and is\n$-\\frac{\\pi}{100}$ for the others. We also explore the effect of\napplying uniform noise to this phase difference.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n_noise = 41 # amount of noise amplitude samples in [0, 4]\ndata = [[]]\n\n# Generate reference\nfor _ in range(n_e):\n data[0].append(np.sin(2 * np.pi * f * t))\n\nA_list = np.linspace(0, 4, n_noise)\n\nfor A in A_list:\n sig = []\n # Generate other signal\n for _ in range(int(n_e / 2)): # phase difference -pi/100\n sig.append(\n np.sin(2 * np.pi * f * t + np.pi / 100 + A * np.random.uniform(-1, 1))\n )\n for _ in range(int(n_e / 2), n_e): # phase difference pi/2\n sig.append(np.sin(2 * np.pi * f * t - np.pi / 2 + A * np.random.uniform(-1, 1)))\n data.append(sig)\n\ndata = np.swapaxes(np.array(data), 0, 1)\n\n# Visualize the data\nfig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True)\nax[0].plot(t[:10], data[0, 0, :10], label=\"Reference\")\nax[0].plot(t[:10], data[1, 1, :10])\n\nax[0].set_title(r\"Phase Lagging ($-\\pi/100$ Phase Difference)\")\nax[0].set_xlabel(\"Time (s)\")\nax[0].set_ylabel(\"Signal\")\nax[0].legend()\n\nax[1].plot(t[:fs], data[0, 0, :fs], label=\"Reference\")\nax[1].plot(t[:fs], data[-1, 1, :fs])\nax[1].set_title(r\"Phase Leading ($\\pi/2$ Phase Difference)\")\nax[1].set_xlabel(\"Time (s)\")\n\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now compute PLI and wPLI\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"conn = []\nindices = ([0] * n_noise, np.arange(1, n_noise + 1))\nfor method in [\"pli\", \"wpli\"]:\n conn.append(\n spectral_connectivity_epochs(\n data,\n method=method,\n sfreq=fs,\n indices=indices,\n fmin=9,\n fmax=11,\n faverage=True,\n ).get_data()[:, 0]\n )\nconn = np.array(conn)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results from the simulation are shown in the figure below. In the case\nwithout noise, the difference between wPLI and PLI is made obvious. In PLI,\nno connectivity is detected, as the $-\\frac{\\pi}{100}$ phase\ndifferences are weighted in the exact same way as the $\\frac{\\pi}{2}$\nrelationships. wPLI is able to avoid the cancellation of the\n$\\frac{\\pi}{2}$ relationships.\n\nAs noise gets added, PLI increases since the $-\\frac{\\pi}{100}$\nrelationships are made positive more often than the $\\frac{\\pi}{2}$\nrelationships are made negative. However, wPLI maintains an advantage in\nits ability to distinguish the underlying structure. Beyond a certain point,\nthe noise dominates any pre-defined structure, and both methods behave\nsimilarly, tending toward $0$. For a more detailed analysis of this\nresult and the properties of wPLI, please refer to (Vinck et al, 2011)\n:footcite:`VinckEtAl2011`.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure()\nplt.plot(A_list, conn[0], \"o-\", label=\"PLI\")\nplt.plot(A_list, conn[1], \"o-\", label=\"wPLI\")\nplt.legend()\nplt.xlabel(\"Noise Amplitude\")\nplt.ylabel(\"Connectivity Measure\")\nplt.title(\"wPLI and PLI Under Increasing Noise\")\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Demo On MEG Data\n\nTo finish this example, we also quickly demonstrate these methods on some\nsample MEG data recorded during visual stimulation.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"data_path = sample.data_path()\nraw_fname = data_path / \"MEG/sample/sample_audvis_filt-0-40_raw.fif\"\nevent_fname = data_path / \"MEG/sample/sample_audvis_filt-0-40_raw-eve.fif\"\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)\n\n\n# Select gradiometers\npicks = mne.pick_types(\n raw.info, meg=\"grad\", eeg=False, stim=False, eog=True, exclude=\"bads\"\n)\n\n# Create epochs\nevent_id, tmin, tmax = 3, -0.2, 1.5 # need a long enough epoch for 5 cycles\nepochs = mne.Epochs(\n raw,\n events,\n event_id,\n tmin,\n tmax,\n picks=picks,\n baseline=(None, 0),\n reject=dict(grad=4000e-13, eog=150e-6),\n)\nepochs.load_data().pick_types(meg=\"grad\") # just keep MEG and no EOG now\n\nfmin, fmax = 4.0, 9.0 # compute connectivity within 4-9 Hz\nsfreq = raw.info[\"sfreq\"] # the sampling frequency\ntmin = 0.0 # exclude the baseline period\n\n# Compute PLI, wPLI, and dPLI\ncon_pli = spectral_connectivity_epochs(\n epochs,\n method=\"pli\",\n mode=\"multitaper\",\n sfreq=sfreq,\n fmin=fmin,\n fmax=fmax,\n faverage=True,\n tmin=tmin,\n mt_adaptive=False,\n n_jobs=1,\n)\n\ncon_wpli = spectral_connectivity_epochs(\n epochs,\n method=\"wpli\",\n mode=\"multitaper\",\n sfreq=sfreq,\n fmin=fmin,\n fmax=fmax,\n faverage=True,\n tmin=tmin,\n mt_adaptive=False,\n n_jobs=1,\n)\n\ncon_dpli = spectral_connectivity_epochs(\n epochs,\n method=\"dpli\",\n mode=\"multitaper\",\n sfreq=sfreq,\n fmin=fmin,\n fmax=fmax,\n faverage=True,\n tmin=tmin,\n mt_adaptive=False,\n n_jobs=1,\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, there is strong connectivity between sensors 190-200 and\nsensors 110-160.\n\nMoreover, after observing the presence of connectivity, dPLI can be used to\nascertain the direction of the phase relationship. Here, it appears that the\ndPLI connectivity in this area is less than $0.5$, and thus sensors\n190-200 are lagging sensors 110-160.\n\nIn keeping with the previous simulation, we can see that wPLI identifies\nstronger connectivity relationships than PLI. This is due to its robustness\nagainst volume conduction effects decreasing the detected connectivity\nstrength, as was mentioned earlier.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axs = plt.subplots(1, 3, figsize=(14, 5), sharey=True)\naxs[0].imshow(con_pli.get_data(\"dense\"), vmin=0, vmax=1)\naxs[0].set_title(\"PLI\")\naxs[0].set_ylabel(\"Sensor 1\")\naxs[0].set_xlabel(\"Sensor 2\")\n\naxs[1].imshow(con_wpli.get_data(\"dense\"), vmin=0, vmax=1)\naxs[1].set_title(\"wPLI\")\naxs[1].set_xlabel(\"Sensor 2\")\n\nim = axs[2].imshow(con_dpli.get_data(\"dense\"), vmin=0, vmax=1)\naxs[2].set_title(\"dPLI\")\naxs[2].set_xlabel(\"Sensor 2\")\n\nfig.colorbar(im, ax=axs.ravel())\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusions\n\nBoth wPLI and dPLI are extensions upon the original PLI method, and provide\ncomplementary information about underlying connectivity.\n\n* To identify the presence of an underlying phase relationship, wPLI is the\n method of choice for most researchers as it provides an improvement in\n robustness over the original PLI method\n* To know the directionality of the connectivity identified by wPLI, dPLI\n should be used\n\nUltimately, these methods work great together, providing a comprehensive\nestimate of phase-based connectivity.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n.. footbibliography::\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}