{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Compute directionality of connectivity with multivariate Granger causality\n\nThis example demonstrates how Granger causality based on state-space models\n:footcite:`BarnettSeth2015` can be used to compute directed connectivity\nbetween sensors in a multivariate manner. Furthermore, the use of time-reversal\nfor improving the robustness of directed connectivity estimates to noise in the\ndata is discussed :footcite:`WinklerEtAl2016`.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Author: Thomas S. Binns \n# License: BSD (3-clause)\n# sphinx_gallery_thumbnail_number = 3"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import mne\nimport numpy as np\nfrom matplotlib import pyplot as plt\nfrom mne.datasets.fieldtrip_cmc import data_path\n\nfrom mne_connectivity import spectral_connectivity_epochs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Background\n\nMultivariate forms of signal analysis allow you to simultaneously consider\nthe activity of multiple signals. In the case of connectivity, the\ninteraction between multiple sensors can be analysed at once, producing a\nsingle connectivity spectrum. This approach brings not only practical\nbenefits (e.g. easier interpretability of results from the dimensionality\nreduction), but can also offer methodological improvements (e.g. enhanced\nsignal-to-noise ratio and reduced bias).\n\nAdditionally, it can be of interest to examine the directionality of\nconnectivity between signals, providing additional clarity to how information\nflows in a system. One such directed measure of connectivity is Granger\ncausality (GC). A signal, $\\boldsymbol{x}$, is said to Granger-cause\nanother signal, $\\boldsymbol{y}$, if information from the past of\n$\\boldsymbol{x}$ improves the prediction of the present of\n$\\boldsymbol{y}$ over the case where only information from the past of\n$\\boldsymbol{y}$ is used. Note: GC does not make any assertions about\nthe true causality between signals.\n\nThe degree to which $\\boldsymbol{x}$ and $\\boldsymbol{y}$ can be\nused to predict one another in a linear model can be quantified using vector\nautoregressive (VAR) models. Considering the simpler case of time domain\nconnectivity, the VAR models are as follows:\n\n$y_t = \\sum_{k=1}^{K} a_k y_{t-k} + \\xi_t^y$ ,\n$Var(\\xi_t^y) := \\Sigma_y$ ,\n\nand $\\boldsymbol{z}_t = \\sum_{k=1}^K \\boldsymbol{A}_k\n\\boldsymbol{z}_{t-k} + \\boldsymbol{\\epsilon}_t$ ,\n$\\boldsymbol{\\Sigma} := \\langle \\boldsymbol{\\epsilon}_t\n\\boldsymbol{\\epsilon}_t^T \\rangle = \\begin{bmatrix} \\Sigma_{xx} & \\Sigma_{xy}\n\\\\ \\Sigma_{yx} & \\Sigma_{yy} \\end{bmatrix}$ ,\n\nrepresenting the reduced and full VAR models, respectively, where: $K$\nis the order of the VAR model, determining the number of lags, $k$,\nused; $\\boldsymbol{Z} := \\begin{bmatrix} \\boldsymbol{x} \\\\\n\\boldsymbol{y} \\end{bmatrix}$; $\\boldsymbol{A}$ is a matrix of\ncoefficients explaining the contribution of past entries of\n$\\boldsymbol{Z}$ to its current value; and $\\xi$ and\n$\\boldsymbol{\\epsilon}$ are the residuals of the VAR models. In this\nway, the information of the signals at time $t$ can be represented as a\nweighted form of the information from the previous timepoints, plus some\nresidual information not encoded in the signals' past. In practice, VAR model\nparameters are computed from an autocovariance sequence generated from the\ntime-series data using the Yule-Walker equations :footcite:`Whittle1963`.\n\nThe residuals, or errors, represent how much information about the present\nstate of the signals is not explained by their past. We can therefore\nestimate how much $\\boldsymbol{x}$ Granger-causes\n$\\boldsymbol{y}$ by comparing the variance of the residuals of the\nreduced VAR model ($\\Sigma_y$; i.e. how much the present of\n$\\boldsymbol{y}$ is not explained by its own past) and of the full VAR\nmodel ($\\Sigma_{yy}$; i.e. how much the present of\n$\\boldsymbol{y}$ is not explained by both its own past and that of\n$\\boldsymbol{x}$):\n\n$F_{x \\rightarrow y} = ln \\Large{(\\frac{\\Sigma_y}{\\Sigma_{yy}})}$ ,\n\nwhere $F$ is the Granger score. For example, if $\\boldsymbol{x}$\ncontains no information about $\\boldsymbol{y}$, the residuals of the\nreduced and full VAR models will be identical, and\n$F_{x \\rightarrow y}$ will naturally be 0, indicating that\ninformation from $\\boldsymbol{x}$ does not flow to\n$\\boldsymbol{y}$. In contrast, if $\\boldsymbol{x}$ does help to\npredict $\\boldsymbol{y}$, the residual of the full model will be\nsmaller than that of the reduced model. $\\Large{\\frac{\\Sigma_y}\n{\\Sigma_{yy}}}$ will therefore be greater than 1, leading to a Granger score\n> 0. Granger scores are bound between $[0, \\infty)$.\n\nThese same principles apply to spectral GC, which provides information about\nthe directionality of connectivity for individual frequencies. For spectral\nGC, the autocovariance sequence is generated from an inverse Fourier\ntransform applied to the cross-spectral density of the signals. Additionally,\na spectral transfer function is used to translate information from the VAR\nmodels back into the frequency domain before computing the final Granger\nscores.\n\nBarnett and Seth (2015) :footcite:`BarnettSeth2015` have defined a\nmultivariate form of spectral GC based on state-space models, enabling the\nestimation of information flow between whole sets of signals simultaneously:\n\n$F_{A \\rightarrow B}(f) = \\Re ln \\Large{(\\frac{\ndet(\\boldsymbol{S}_{BB}(f))}{det(\\boldsymbol{S}_{BB}(f) -\n\\boldsymbol{H}_{BA}(f) \\boldsymbol{\\Sigma}_{AA \\lvert B}\n\\boldsymbol{H}_{BA}^*(f))})}$ ,\n\nwhere: $A$ and $B$ are the seeds and targets, respectively;\n$f$ is a given frequency; $\\boldsymbol{H}$ is the spectral\ntransfer function; $\\boldsymbol{\\Sigma}$ is the innovations form\nresiduals' covariance matrix of the state-space model; $\\boldsymbol{S}$\nis $\\boldsymbol{\\Sigma}$ transformed by $\\boldsymbol{H}$; and\n$\\boldsymbol{\\Sigma}_{IJ \\lvert K} := \\boldsymbol{\\Sigma}_{IJ} -\n\\boldsymbol{\\Sigma}_{IK} \\boldsymbol{\\Sigma}_{KK}^{-1}\n\\boldsymbol{\\Sigma}_{KJ}$, representing a partial covariance matrix. The same\nprinciples apply as before: a numerator greater than the denominator means\nthat information from the seed signals aids the prediction of activity in the\ntarget signals, leading to a Granger score > 0.\n\nThere are several benefits to a state-space approach for computing GC:\ncompared to traditional autoregressive-based approaches, the use of\nstate-space models offers reduced statistical bias and increased statistical\npower; furthermore, the dimensionality reduction offered by the multivariate\nnature of the approach can aid in the interpretability and subsequent\nanalysis of the results.\n\nTo demonstrate the use of GC for estimating directed connectivity, we start\nby loading some example MEG data and dividing it into two-second-long epochs.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"raw = mne.io.read_raw_ctf(data_path() / \"SubjectCMC.ds\")\nraw.pick(\"mag\")\nraw.crop(50.0, 110.0).load_data()\nraw.notch_filter(50)\nraw.resample(100)\n\nepochs = mne.make_fixed_length_epochs(raw, duration=2.0).load_data()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will focus on connectivity between sensors over the parietal and occipital\ncortices, with 20 parietal sensors designated as group A, and 22 occipital\nsensors designated as group B.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# parietal sensors\nsignals_a = [\n idx\n for idx, ch_info in enumerate(epochs.info[\"chs\"])\n if ch_info[\"ch_name\"][2] == \"P\"\n]\n# occipital sensors\nsignals_b = [\n idx\n for idx, ch_info in enumerate(epochs.info[\"chs\"])\n if ch_info[\"ch_name\"][2] == \"O\"\n]\n\nindices_ab = (np.array([signals_a]), np.array([signals_b])) # A => B\nindices_ba = (np.array([signals_b]), np.array([signals_a])) # B => A\n\n# compute Granger causality\ngc_ab = spectral_connectivity_epochs(\n epochs,\n method=[\"gc\"],\n indices=indices_ab,\n fmin=5,\n fmax=30,\n rank=(np.array([5]), np.array([5])),\n gc_n_lags=20,\n) # A => B\ngc_ba = spectral_connectivity_epochs(\n epochs,\n method=[\"gc\"],\n indices=indices_ba,\n fmin=5,\n fmax=30,\n rank=(np.array([5]), np.array([5])),\n gc_n_lags=20,\n) # B => A\nfreqs = gc_ab.freqs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting the results, we see that there is a flow of information from our\nparietal sensors (group A) to our occipital sensors (group B) with a\nnoticeable peak at ~8 Hz, and smaller peaks at 18 and 26 Hz.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axis = plt.subplots(1, 1)\naxis.plot(freqs, gc_ab.get_data()[0], linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Connectivity (A.U.)\")\nfig.suptitle(\"GC: [A => B]\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Drivers and receivers: analysing the net direction of information flow\n\nAlthough analysing connectivity in a given direction can be of interest,\nthere may exist a bidirectional relationship between signals. In such cases,\nidentifying the signals that dominate information flow (the drivers) may be\ndesired. For this, we can simply subtract the Granger scores in the opposite\ndirection, giving us the net GC score:\n\n$F_{A \\rightarrow B}^{net} := F_{A \\rightarrow B} -\nF_{B \\rightarrow A}$.\n\nDoing so, we see that the flow of information across the spectrum remains\ndominant from parietal to occipital sensors (indicated by the positive-valued\nGranger scores), with similar peaks around 10, 18, and 26 Hz.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"net_gc = gc_ab.get_data() - gc_ba.get_data() # [A => B] - [B => A]\n\nfig, axis = plt.subplots(1, 1)\naxis.plot((freqs[0], freqs[-1]), (0, 0), linewidth=2, linestyle=\"--\", color=\"k\")\naxis.plot(freqs, net_gc[0], linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Connectivity (A.U.)\")\nfig.suptitle(\"Net GC: [A => B] - [B => A]\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Improving the robustness of connectivity estimates with time-reversal\n\nOne limitation of GC methods is the risk of connectivity estimates being\ncontaminated with noise. For instance, consider the case where, due to\nvolume conduction, multiple sensors detect activity from the same source.\nNaturally, information recorded at these sensors mutually help to predict\nthe activity of one another, leading to spurious estimates of directed\nconnectivity which one may incorrectly attribute to information flow between\ndifferent brain regions. On the other hand, even if there is no source\nmixing, the presence of correlated noise between sensors can similarly bias\ndirected connectivity estimates.\n\nTo address this issue, Haufe *et al.* (2013) :footcite:`HaufeEtAl2013`\npropose contrasting causality scores obtained on the original time-series to\nthose obtained on the reversed time-series. The idea behind this approach is\nas follows: if temporal order is crucial in distinguishing a driver from a\nrecipient, then reversing the temporal order should reduce, if not flip, an\nestimate of directed connectivity. In practice, time-reversal is implemented\nas a transposition of the autocovariance sequence used to compute GC\n:footcite:`HaufeEtAl2012`.\n\nSeveral studies have shown that that such an approach can reduce the degree\nof false-positive connectivity estimates (even performing favourably against\nother methods such as the phase slope index) :footcite:`VinckEtAl2015` and\nretain the ability to correctly identify the net direction of information\nflow akin to net GC :footcite:`WinklerEtAl2016,HaufeEtAl2013`. This approach\nis termed time-reversed GC (TRGC):\n\n$\\tilde{D}_{A \\rightarrow B}^{net} := F_{A \\rightarrow B}^{net} -\nF_{\\tilde{A} \\rightarrow \\tilde{B}}^{net}$ ,\n\nwhere $\\sim$ represents time-reversal, and:\n\n$F_{\\tilde{A} \\rightarrow \\tilde{B}}^{net} := F_{\\tilde{A} \\rightarrow\n\\tilde{B}} - F_{\\tilde{B} \\rightarrow \\tilde{A}}$.\n\nGC on time-reversed signals can be computed simply with ``method=['gc_tr']``,\nwhich will perform the time-reversal of the signals for the end-user. Note\nthat **time-reversed results should only be interpreted in the context of net\nresults**, i.e. with $\\tilde{D}_{A \\rightarrow B}^{net}$. In the\nexample below, notice how the outputs are not used directly, but rather used\nto produce net scores of the time-reversed signals. The net scores of the\ntime-reversed signals can then be subtracted from the net scores of the\noriginal signals to produce the final TRGC scores.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# compute GC on time-reversed signals\ngc_tr_ab = spectral_connectivity_epochs(\n epochs,\n method=[\"gc_tr\"],\n indices=indices_ab,\n fmin=5,\n fmax=30,\n rank=(np.array([5]), np.array([5])),\n gc_n_lags=20,\n) # TR[A => B]\ngc_tr_ba = spectral_connectivity_epochs(\n epochs,\n method=[\"gc_tr\"],\n indices=indices_ba,\n fmin=5,\n fmax=30,\n rank=(np.array([5]), np.array([5])),\n gc_n_lags=20,\n) # TR[B => A]\n\n# compute net GC on time-reversed signals (TR[A => B] - TR[B => A])\nnet_gc_tr = gc_tr_ab.get_data() - gc_tr_ba.get_data()\n\n# compute TRGC\ntrgc = net_gc - net_gc_tr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting the TRGC results reveals a very different picture compared to net\nGC. For one, there is now a dominance of information flow ~6 Hz from\noccipital to parietal sensors (indicated by the negative-valued Granger\nscores). Additionally, the peak ~10 Hz is less dominant in the spectrum, with\nparietal to occipital information flow between 13-20 Hz being much more\nprominent. The stark difference between net GC and TRGC results indicates\nthat the net GC spectrum was contaminated by spurious connectivity resulting\nfrom source mixing or correlated noise in the recordings. Altogether, the use\nof TRGC instead of net GC is generally advised.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axis = plt.subplots(1, 1)\naxis.plot((freqs[0], freqs[-1]), (0, 0), linewidth=2, linestyle=\"--\", color=\"k\")\naxis.plot(freqs, trgc[0], linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Connectivity (A.U.)\")\nfig.suptitle(\"TRGC: net[A => B] - net time-reversed[A => B]\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Controlling spectral smoothing with the number of lags\n\nOne important parameter when computing GC is the number of lags used when\ncomputing the VAR model. A lower number of lags reduces the computational\ncost, but in the context of spectral GC, leads to a smoothing of Granger\nscores across frequencies. The number of lags can be specified using the\n``gc_n_lags`` parameter. The default value is 40, however there is no correct\nnumber of lags to use when computing GC. Instead, you have to use your own\nbest judgement of whether or not your Granger scores look overly smooth.\n\nBelow is a comparison of Granger scores computed with a different number of\nlags. In the above examples we used 20 lags, which we will compare to Granger\nscores computed with 60 lags. As you can see, the spectra of Granger scores\ncomputed with 60 lags is noticeably less smooth, but it does share the same\noverall pattern.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"gc_ab_60 = spectral_connectivity_epochs(\n epochs,\n method=[\"gc\"],\n indices=indices_ab,\n fmin=5,\n fmax=30,\n rank=(np.array([5]), np.array([5])),\n gc_n_lags=60,\n) # A => B\n\nfig, axis = plt.subplots(1, 1)\naxis.plot(freqs, gc_ab.get_data()[0], linewidth=2, label=\"20 lags\")\naxis.plot(freqs, gc_ab_60.get_data()[0], linewidth=2, label=\"60 lags\")\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Connectivity (A.U.)\")\naxis.legend()\nfig.suptitle(\"GC: [A => B]\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Handling high-dimensional data\n\nAn important issue to consider when computing multivariate GC is that the\ndata GC is computed on should not be rank deficient (i.e. must have full\nrank). More specifically, the autocovariance matrix must not be singular or\nclose to singular.\n\nIn the case that your data is not full rank and ``rank`` is left as ``None``,\nan automatic rank computation is performed and an appropriate degree of\ndimensionality reduction will be enforced. The rank of the data is determined\nby computing the singular values of the data and finding those within a\nfactor of $1e^{-6}$ relative to the largest singular value.\n\nWhilst unlikely, there may be scenarios in which this threshold is too\nlenient. In these cases, you should inspect the singular values of your data\nto identify an appropriate degree of dimensionality reduction to perform,\nwhich you can then specify manually using the ``rank`` argument. The code\nbelow shows one possible approach for finding an appropriate rank of\nclose-to-singular data with a more conservative threshold.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# gets the singular values of the data\ns = np.linalg.svd(raw.get_data(), compute_uv=False)\n# finds how many singular values are 'close' to the largest singular value\nrank = np.count_nonzero(s >= s[0] * 1e-4) # 1e-4 is the 'closeness' criteria"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nonethless, even in situations where you specify an appropriate rank, it is\nnot guaranteed that the subsequently-computed autocovariance sequence will\nretain this non-singularity (this can depend on, e.g. the number of lags).\nHence, you may also encounter situations where you have to specify a rank\nless than that of your data to ensure that the autocovariance sequence is\nnon-singular.\n\nIn the above examples, notice how a rank of 5 was given, despite there being\n20 channels in the seeds and targets. Attempting to compute GC on the\noriginal data would not succeed, given that the resulting autocovariance\nsequence is singular, as the example below shows.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"try:\n spectral_connectivity_epochs(\n epochs,\n method=[\"gc\"],\n indices=indices_ab,\n fmin=5,\n fmax=30,\n rank=None,\n gc_n_lags=20,\n verbose=False,\n ) # A => B\n print(\"Success!\")\nexcept RuntimeError as error:\n print(\"\\nCaught the following error:\\n\" + repr(error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Rigorous checks are implemented to identify any such instances which would\notherwise cause the GC computation to produce erroneous results. You can\ntherefore be confident as an end-user that these cases will be caught.\n\nFinally, when comparing GC scores across recordings, **it is highly\nrecommended to estimate connectivity from the same number of channels (or\nequally from the same degree of rank subspace projection)** to avoid biases\nin connectivity estimates. Bias can be avoided by specifying a consistent\nrank subspace to project to using the ``rank`` argument, standardising your\nconnectivity estimates regardless of changes in e.g. the number of channels\nacross recordings. Note that this does not refer to the number of seeds and\ntargets *within* a connection being identical, rather to the number of seeds\nand targets *across* connections.\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
}