Compute directionality of connectivity with multivariate Granger causality

This example demonstrates how Granger causality based on state-space models [1] can be used to compute directed connectivity between sensors in a multivariate manner. Furthermore, the use of time-reversal for improving the robustness of directed connectivity estimates to noise in the data is discussed [2].

# Author: Thomas S. Binns <t.s.binns@outlook.com>
# License: BSD (3-clause)
# sphinx_gallery_failing_thumbnail = False
# sphinx_gallery_thumbnail_number = 3

import mne
import numpy as np
from matplotlib import pyplot as plt
from mne.datasets.fieldtrip_cmc import data_path

from mne_connectivity import spectral_connectivity_epochs

Background

Multivariate forms of signal analysis allow you to simultaneously consider the activity of multiple signals. In the case of connectivity, the interaction between multiple sensors can be analysed at once, producing a single connectivity spectrum. This approach brings not only practical benefits (e.g. easier interpretability of results from the dimensionality reduction), but can also offer methodological improvements (e.g. enhanced signal-to-noise ratio and reduced bias).

Additionally, it can be of interest to examine the directionality of connectivity between signals, providing additional clarity to how information flows in a system. One such directed measure of connectivity is Granger causality (GC). A signal, \(\boldsymbol{x}\), is said to Granger-cause another signal, \(\boldsymbol{y}\), if information from the past of \(\boldsymbol{x}\) improves the prediction of the present of \(\boldsymbol{y}\) over the case where only information from the past of \(\boldsymbol{y}\) is used. Note: GC does not make any assertions about the true causality between signals.

The degree to which \(\boldsymbol{x}\) and \(\boldsymbol{y}\) can be used to predict one another in a linear model can be quantified using vector autoregressive (VAR) models. Considering the simpler case of time domain connectivity, the VAR models are as follows:

\(y_t = \sum_{k=1}^{K} a_k y_{t-k} + \xi_t^y\) , \(Var(\xi_t^y) := \Sigma_y\) ,

and \(\boldsymbol{z}_t = \sum_{k=1}^K \boldsymbol{A}_k \boldsymbol{z}_{t-k} + \boldsymbol{\epsilon}_t\) , \(\boldsymbol{\Sigma} := \langle \boldsymbol{\epsilon}_t \boldsymbol{\epsilon}_t^T \rangle = \begin{bmatrix} \Sigma_{xx} & \Sigma_{xy} \\ \Sigma_{yx} & \Sigma_{yy} \end{bmatrix}\) ,

representing the reduced and full VAR models, respectively, where: \(K\) is the order of the VAR model, determining the number of lags, \(k\), used; \(\boldsymbol{Z} := \begin{bmatrix} \boldsymbol{x} \\ \boldsymbol{y} \end{bmatrix}\); \(\boldsymbol{A}\) is a matrix of coefficients explaining the contribution of past entries of \(\boldsymbol{Z}\) to its current value; and \(\xi\) and \(\boldsymbol{\epsilon}\) are the residuals of the VAR models. In this way, the information of the signals at time \(t\) can be represented as a weighted form of the information from the previous timepoints, plus some residual information not encoded in the signals’ past. In practice, VAR model parameters are computed from an autocovariance sequence generated from the time-series data using the Yule-Walker equations [3].

The residuals, or errors, represent how much information about the present state of the signals is not explained by their past. We can therefore estimate how much \(\boldsymbol{x}\) Granger-causes \(\boldsymbol{y}\) by comparing the variance of the residuals of the reduced VAR model (\(\Sigma_y\); i.e. how much the present of \(\boldsymbol{y}\) is not explained by its own past) and of the full VAR model (\(\Sigma_{yy}\); i.e. how much the present of \(\boldsymbol{y}\) is not explained by both its own past and that of \(\boldsymbol{x}\)):

\(F_{x \rightarrow y} = ln \Large{(\frac{\Sigma_y}{\Sigma_{yy}})}\) ,

where \(F\) is the Granger score. For example, if \(\boldsymbol{x}\) contains no information about \(\boldsymbol{y}\), the residuals of the reduced and full VAR models will be identical, and \(F_{x \rightarrow y}\) will naturally be 0, indicating that information from \(\boldsymbol{x}\) does not flow to \(\boldsymbol{y}\). In contrast, if \(\boldsymbol{x}\) does help to predict \(\boldsymbol{y}\), the residual of the full model will be smaller than that of the reduced model. \(\Large{\frac{\Sigma_y} {\Sigma_{yy}}}\) will therefore be greater than 1, leading to a Granger score > 0. Granger scores are bound between \([0, \infty)\).

These same principles apply to spectral GC, which provides information about the directionality of connectivity for individual frequencies. For spectral GC, the autocovariance sequence is generated from an inverse Fourier transform applied to the cross-spectral density of the signals. Additionally, a spectral transfer function is used to translate information from the VAR models back into the frequency domain before computing the final Granger scores.

Barnett and Seth (2015) [1] have defined a multivariate form of spectral GC based on state-space models, enabling the estimation of information flow between whole sets of signals simultaneously:

\(F_{A \rightarrow B}(f) = \Re ln \Large{(\frac{ det(\boldsymbol{S}_{BB}(f))}{det(\boldsymbol{S}_{BB}(f) - \boldsymbol{H}_{BA}(f) \boldsymbol{\Sigma}_{AA \lvert B} \boldsymbol{H}_{BA}^*(f))})}\) ,

where: \(A\) and \(B\) are the seeds and targets, respectively; \(f\) is a given frequency; \(\boldsymbol{H}\) is the spectral transfer function; \(\boldsymbol{\Sigma}\) is the innovations form residuals’ covariance matrix of the state-space model; \(\boldsymbol{S}\) is \(\boldsymbol{\Sigma}\) transformed by \(\boldsymbol{H}\); and \(\boldsymbol{\Sigma}_{IJ \lvert K} := \boldsymbol{\Sigma}_{IJ} - \boldsymbol{\Sigma}_{IK} \boldsymbol{\Sigma}_{KK}^{-1} \boldsymbol{\Sigma}_{KJ}\), representing a partial covariance matrix. The same principles apply as before: a numerator greater than the denominator means that information from the seed signals aids the prediction of activity in the target signals, leading to a Granger score > 0.

There are several benefits to a state-space approach for computing GC: compared to traditional autoregressive-based approaches, the use of state-space models offers reduced statistical bias and increased statistical power; furthermore, the dimensionality reduction offered by the multivariate nature of the approach can aid in the interpretability and subsequent analysis of the results.

To demonstrate the use of GC for estimating directed connectivity, we start by loading some example MEG data and dividing it into two-second-long epochs.

raw = mne.io.read_raw_ctf(data_path() / "SubjectCMC.ds")
raw.pick("mag")
raw.crop(50.0, 110.0).load_data()
raw.notch_filter(50)
raw.resample(100)

epochs = mne.make_fixed_length_epochs(raw, duration=2.0).load_data()

  0%|                                               | 0.00/329M [00:00<?, ?B/s]
  1%|▌                                     | 4.83M/329M [00:00<00:06, 48.2MB/s]
  3%|█▏                                    | 10.6M/329M [00:00<00:05, 53.9MB/s]
  5%|█▉                                    | 16.7M/329M [00:00<00:05, 57.1MB/s]
  7%|██▋                                   | 22.9M/329M [00:00<00:05, 59.0MB/s]
  9%|███▎                                  | 29.2M/329M [00:00<00:04, 60.5MB/s]
 11%|████▏                                 | 36.5M/329M [00:00<00:04, 64.5MB/s]
 13%|████▉                                 | 43.1M/329M [00:00<00:04, 65.1MB/s]
 15%|█████▊                                | 50.5M/329M [00:00<00:04, 67.8MB/s]
 18%|██████▊                               | 58.9M/329M [00:00<00:03, 73.1MB/s]
 20%|███████▋                              | 67.0M/329M [00:01<00:03, 75.6MB/s]
 23%|████████▋                             | 75.1M/329M [00:01<00:03, 77.1MB/s]
 25%|█████████▋                            | 83.4M/329M [00:01<00:03, 78.9MB/s]
 28%|██████████▌                           | 91.8M/329M [00:01<00:02, 80.5MB/s]
 30%|███████████▌                          | 99.9M/329M [00:01<00:02, 79.2MB/s]
 33%|████████████▊                          | 108M/329M [00:01<00:02, 80.7MB/s]
 35%|█████████████▊                         | 116M/329M [00:01<00:02, 80.5MB/s]
 38%|██████████████▋                        | 124M/329M [00:01<00:02, 80.5MB/s]
 40%|███████████████▋                       | 132M/329M [00:01<00:02, 79.9MB/s]
 43%|████████████████▋                      | 140M/329M [00:01<00:02, 79.4MB/s]
 45%|█████████████████▋                     | 149M/329M [00:02<00:02, 80.6MB/s]
 48%|██████████████████▌                    | 157M/329M [00:02<00:02, 78.5MB/s]
 50%|███████████████████▌                   | 165M/329M [00:02<00:02, 75.7MB/s]
 52%|████████████████████▍                  | 172M/329M [00:02<00:02, 71.1MB/s]
 55%|█████████████████████▎                 | 179M/329M [00:02<00:02, 67.3MB/s]
 57%|██████████████████████                 | 186M/329M [00:02<00:02, 66.1MB/s]
 59%|██████████████████████▊                | 193M/329M [00:02<00:02, 66.4MB/s]
 61%|███████████████████████▋               | 200M/329M [00:02<00:02, 62.8MB/s]
 63%|████████████████████████▍              | 206M/329M [00:02<00:02, 57.5MB/s]
 64%|█████████████████████████              | 212M/329M [00:03<00:02, 54.5MB/s]
 66%|█████████████████████████▋             | 217M/329M [00:03<00:02, 52.5MB/s]
 68%|██████████████████████████▍            | 223M/329M [00:03<00:02, 51.0MB/s]
 69%|██████████████████████████▉            | 228M/329M [00:03<00:02, 50.1MB/s]
 71%|███████████████████████████▌           | 233M/329M [00:03<00:01, 49.3MB/s]
 72%|████████████████████████████▏          | 238M/329M [00:03<00:01, 48.9MB/s]
 74%|████████████████████████████▋          | 243M/329M [00:03<00:01, 48.7MB/s]
 76%|█████████████████████████████▍         | 249M/329M [00:03<00:01, 53.0MB/s]
 78%|██████████████████████████████▎        | 255M/329M [00:03<00:01, 56.3MB/s]
 79%|██████████████████████████████▉        | 262M/329M [00:04<00:01, 57.7MB/s]
 81%|███████████████████████████████▋       | 268M/329M [00:04<00:01, 59.5MB/s]
 83%|████████████████████████████████▌      | 275M/329M [00:04<00:00, 62.5MB/s]
 85%|█████████████████████████████████▎     | 281M/329M [00:04<00:00, 61.7MB/s]
 87%|██████████████████████████████████     | 287M/329M [00:04<00:00, 59.9MB/s]
 89%|██████████████████████████████████▊    | 293M/329M [00:04<00:00, 58.1MB/s]
 91%|███████████████████████████████████▍   | 299M/329M [00:04<00:00, 56.6MB/s]
 93%|████████████████████████████████████   | 305M/329M [00:04<00:00, 56.1MB/s]
 94%|████████████████████████████████████▊  | 310M/329M [00:04<00:00, 54.3MB/s]
 96%|█████████████████████████████████████▍ | 316M/329M [00:04<00:00, 55.1MB/s]
 98%|██████████████████████████████████████ | 322M/329M [00:05<00:00, 54.6MB/s]
 99%|██████████████████████████████████████▊| 327M/329M [00:05<00:00, 54.9MB/s]
  0%|                                               | 0.00/329M [00:00<?, ?B/s]
100%|████████████████████████████████████████| 329M/329M [00:00<00:00, 900GB/s]
Download complete in 13s (314.0 MB)
ds directory : /home/circleci/mne_data/MNE-fieldtrip_cmc-data/SubjectCMC.ds
    res4 data read.
    hc data read.
    Separate EEG position data file not present.
    Quaternion matching (desired vs. transformed):
       0.33   78.32    0.00 mm <->    0.33   78.32    0.00 mm (orig :  -71.62   40.46 -256.48 mm) diff =    0.000 mm
      -0.33  -78.32   -0.00 mm <->   -0.33  -78.32   -0.00 mm (orig :   39.27  -70.16 -258.60 mm) diff =    0.000 mm
     114.65    0.00   -0.00 mm <->  114.65    0.00   -0.00 mm (orig :   64.35   66.64 -262.01 mm) diff =    0.000 mm
    Coordinate transformations established.
    Polhemus data for 3 HPI coils added
    Device coordinate locations for 3 HPI coils added
Picked positions of 4 EEG channels from channel info
    4 EEG locations added to Polhemus data.
    Measurement info composed.
Finding samples for /home/circleci/mne_data/MNE-fieldtrip_cmc-data/SubjectCMC.ds/SubjectCMC.meg4:
    System clock channel is available, checking which samples are valid.
    75 x 12000 = 911610 samples from 191 chs
    390 samples omitted at the end
Current compensation grade : 0
Removing 5 compensators from info because not all compensation channels were picked.
Reading 0 ... 72000  =      0.000 ...    60.000 secs...
Filtering raw data in 1 contiguous segment
Setting up band-stop filter from 49 - 51 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 49.38
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 49.12 Hz)
- Upper passband edge: 50.62 Hz
- Upper transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 50.88 Hz)
- Filter length: 7921 samples (6.601 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.3s
Not setting metadata
30 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 30 events and 200 original time points ...
0 bad epochs dropped

We will focus on connectivity between sensors over the parietal and occipital cortices, with 20 parietal sensors designated as group A, and 22 occipital sensors designated as group B.

# parietal sensors
signals_a = [
    idx
    for idx, ch_info in enumerate(epochs.info["chs"])
    if ch_info["ch_name"][2] == "P"
]
# occipital sensors
signals_b = [
    idx
    for idx, ch_info in enumerate(epochs.info["chs"])
    if ch_info["ch_name"][2] == "O"
]

indices_ab = (np.array([signals_a]), np.array([signals_b]))  # A => B
indices_ba = (np.array([signals_b]), np.array([signals_a]))  # B => A

# compute Granger causality
gc_ab = spectral_connectivity_epochs(
    epochs,
    method=["gc"],
    indices=indices_ab,
    fmin=5,
    fmax=30,
    rank=(np.array([5]), np.array([5])),
    gc_n_lags=20,
)  # A => B
gc_ba = spectral_connectivity_epochs(
    epochs,
    method=["gc"],
    indices=indices_ba,
    fmin=5,
    fmax=30,
    rank=(np.array([5]), np.array([5])),
    gc_n_lags=20,
)  # B => A
freqs = gc_ab.freqs

Adding metadata with 3 columns
Connectivity computation...
    using t=0.000s..1.990s for estimation (200 points)
    frequencies: 5.0Hz..30.0Hz (51 points)
    computing connectivity for 1 connections
    Using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: GC
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
Computing GC for connection 1 of 1

  0%|          | frequency blocks : 0/51 [00:00<?,       ?it/s]
 14%|█▎        | frequency blocks : 7/51 [00:00<00:00,  423.34it/s]
 31%|███▏      | frequency blocks : 16/51 [00:00<00:00,  469.71it/s]
 55%|█████▍    | frequency blocks : 28/51 [00:00<00:00,  539.26it/s]
 75%|███████▍  | frequency blocks : 38/51 [00:00<00:00,  557.30it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00,  619.07it/s]
[Connectivity computation done]
Replacing existing metadata with 3 columns
Connectivity computation...
    using t=0.000s..1.990s for estimation (200 points)
    frequencies: 5.0Hz..30.0Hz (51 points)
    computing connectivity for 1 connections
    Using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: GC
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
Computing GC for connection 1 of 1

  0%|          | frequency blocks : 0/51 [00:00<?,       ?it/s]
 20%|█▉        | frequency blocks : 10/51 [00:00<00:00,  599.22it/s]
 43%|████▎     | frequency blocks : 22/51 [00:00<00:00,  663.28it/s]
 65%|██████▍   | frequency blocks : 33/51 [00:00<00:00,  665.97it/s]
 88%|████████▊ | frequency blocks : 45/51 [00:00<00:00,  683.29it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00,  693.34it/s]
[Connectivity computation done]

Plotting the results, we see that there is a flow of information from our parietal sensors (group A) to our occipital sensors (group B) with a noticeable peak at ~8 Hz, and smaller peaks at 18 and 26 Hz.

fig, axis = plt.subplots(1, 1)
axis.plot(freqs, gc_ab.get_data()[0], linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Connectivity (A.U.)")
fig.suptitle("GC: [A => B]")

Text(0.5, 0.98, 'GC: [A => B]')

Drivers and receivers: analysing the net direction of information flow

Although analysing connectivity in a given direction can be of interest, there may exist a bidirectional relationship between signals. In such cases, identifying the signals that dominate information flow (the drivers) may be desired. For this, we can simply subtract the Granger scores in the opposite direction, giving us the net GC score:

\(F_{A \rightarrow B}^{net} := F_{A \rightarrow B} - F_{B \rightarrow A}\).

Doing so, we see that the flow of information across the spectrum remains dominant from parietal to occipital sensors (indicated by the positive-valued Granger scores), with similar peaks around 10, 18, and 26 Hz.

net_gc = gc_ab.get_data() - gc_ba.get_data()  # [A => B] - [B => A]

fig, axis = plt.subplots(1, 1)
axis.plot((freqs[0], freqs[-1]), (0, 0), linewidth=2, linestyle="--", color="k")
axis.plot(freqs, net_gc[0], linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Connectivity (A.U.)")
fig.suptitle("Net GC: [A => B] - [B => A]")

Text(0.5, 0.98, 'Net GC: [A => B] - [B => A]')

Improving the robustness of connectivity estimates with time-reversal

One limitation of GC methods is the risk of connectivity estimates being contaminated with noise. For instance, consider the case where, due to volume conduction, multiple sensors detect activity from the same source. Naturally, information recorded at these sensors mutually help to predict the activity of one another, leading to spurious estimates of directed connectivity which one may incorrectly attribute to information flow between different brain regions. On the other hand, even if there is no source mixing, the presence of correlated noise between sensors can similarly bias directed connectivity estimates.

To address this issue, Haufe et al. (2013) [4] propose contrasting causality scores obtained on the original time-series to those obtained on the reversed time-series. The idea behind this approach is as follows: if temporal order is crucial in distinguishing a driver from a recipient, then reversing the temporal order should reduce, if not flip, an estimate of directed connectivity. In practice, time-reversal is implemented as a transposition of the autocovariance sequence used to compute GC [5].

Several studies have shown that that such an approach can reduce the degree of false-positive connectivity estimates (even performing favourably against other methods such as the phase slope index) [6] and retain the ability to correctly identify the net direction of information flow akin to net GC [2][4]. This approach is termed time-reversed GC (TRGC):

\(\tilde{D}_{A \rightarrow B}^{net} := F_{A \rightarrow B}^{net} - F_{\tilde{A} \rightarrow \tilde{B}}^{net}\) ,

where \(\sim\) represents time-reversal, and:

\(F_{\tilde{A} \rightarrow \tilde{B}}^{net} := F_{\tilde{A} \rightarrow \tilde{B}} - F_{\tilde{B} \rightarrow \tilde{A}}\).

GC on time-reversed signals can be computed simply with method=['gc_tr'], which will perform the time-reversal of the signals for the end-user. Note that time-reversed results should only be interpreted in the context of net results, i.e. with \(\tilde{D}_{A \rightarrow B}^{net}\). In the example below, notice how the outputs are not used directly, but rather used to produce net scores of the time-reversed signals. The net scores of the time-reversed signals can then be subtracted from the net scores of the original signals to produce the final TRGC scores.

# compute GC on time-reversed signals
gc_tr_ab = spectral_connectivity_epochs(
    epochs,
    method=["gc_tr"],
    indices=indices_ab,
    fmin=5,
    fmax=30,
    rank=(np.array([5]), np.array([5])),
    gc_n_lags=20,
)  # TR[A => B]
gc_tr_ba = spectral_connectivity_epochs(
    epochs,
    method=["gc_tr"],
    indices=indices_ba,
    fmin=5,
    fmax=30,
    rank=(np.array([5]), np.array([5])),
    gc_n_lags=20,
)  # TR[B => A]

# compute net GC on time-reversed signals (TR[A => B] - TR[B => A])
net_gc_tr = gc_tr_ab.get_data() - gc_tr_ba.get_data()

# compute TRGC
trgc = net_gc - net_gc_tr

Replacing existing metadata with 3 columns
Connectivity computation...
    using t=0.000s..1.990s for estimation (200 points)
    frequencies: 5.0Hz..30.0Hz (51 points)
    computing connectivity for 1 connections
    Using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: GC time-reversed
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
Computing GC time-reversed for connection 1 of 1

  0%|          | frequency blocks : 0/51 [00:00<?,       ?it/s]
 24%|██▎       | frequency blocks : 12/51 [00:00<00:00,  689.94it/s]
 49%|████▉     | frequency blocks : 25/51 [00:00<00:00,  749.32it/s]
 69%|██████▊   | frequency blocks : 35/51 [00:00<00:00,  673.43it/s]
 98%|█████████▊| frequency blocks : 50/51 [00:00<00:00,  720.98it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00,  707.34it/s]
[Connectivity computation done]
Replacing existing metadata with 3 columns
Connectivity computation...
    using t=0.000s..1.990s for estimation (200 points)
    frequencies: 5.0Hz..30.0Hz (51 points)
    computing connectivity for 1 connections
    Using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: GC time-reversed
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
Computing GC time-reversed for connection 1 of 1

  0%|          | frequency blocks : 0/51 [00:00<?,       ?it/s]
 18%|█▊        | frequency blocks : 9/51 [00:00<00:00,  546.13it/s]
 41%|████      | frequency blocks : 21/51 [00:00<00:00,  645.64it/s]
 63%|██████▎   | frequency blocks : 32/51 [00:00<00:00,  652.64it/s]
 84%|████████▍ | frequency blocks : 43/51 [00:00<00:00,  649.23it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00,  613.91it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00,  613.32it/s]
[Connectivity computation done]

Plotting the TRGC results reveals a very different picture compared to net GC. For one, there is now a dominance of information flow ~6 Hz from occipital to parietal sensors (indicated by the negative-valued Granger scores). Additionally, the peak ~10 Hz is less dominant in the spectrum, with parietal to occipital information flow between 13-20 Hz being much more prominent. The stark difference between net GC and TRGC results indicates that the net GC spectrum was contaminated by spurious connectivity resulting from source mixing or correlated noise in the recordings. Altogether, the use of TRGC instead of net GC is generally advised.

fig, axis = plt.subplots(1, 1)
axis.plot((freqs[0], freqs[-1]), (0, 0), linewidth=2, linestyle="--", color="k")
axis.plot(freqs, trgc[0], linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Connectivity (A.U.)")
fig.suptitle("TRGC: net[A => B] - net time-reversed[A => B]")

Text(0.5, 0.98, 'TRGC: net[A => B] - net time-reversed[A => B]')

Controlling spectral smoothing with the number of lags

One important parameter when computing GC is the number of lags used when computing the VAR model. A lower number of lags reduces the computational cost, but in the context of spectral GC, leads to a smoothing of Granger scores across frequencies. The number of lags can be specified using the gc_n_lags parameter. The default value is 40, however there is no correct number of lags to use when computing GC. Instead, you have to use your own best judgement of whether or not your Granger scores look overly smooth.

Below is a comparison of Granger scores computed with a different number of lags. In the above examples we used 20 lags, which we will compare to Granger scores computed with 60 lags. As you can see, the spectra of Granger scores computed with 60 lags is noticeably less smooth, but it does share the same overall pattern.

gc_ab_60 = spectral_connectivity_epochs(
    epochs,
    method=["gc"],
    indices=indices_ab,
    fmin=5,
    fmax=30,
    rank=(np.array([5]), np.array([5])),
    gc_n_lags=60,
)  # A => B

fig, axis = plt.subplots(1, 1)
axis.plot(freqs, gc_ab.get_data()[0], linewidth=2, label="20 lags")
axis.plot(freqs, gc_ab_60.get_data()[0], linewidth=2, label="60 lags")
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Connectivity (A.U.)")
axis.legend()
fig.suptitle("GC: [A => B]")

Replacing existing metadata with 3 columns
Connectivity computation...
    using t=0.000s..1.990s for estimation (200 points)
    frequencies: 5.0Hz..30.0Hz (51 points)
    computing connectivity for 1 connections
    Using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: GC
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
Computing GC for connection 1 of 1

  0%|          | frequency blocks : 0/51 [00:00<?,       ?it/s]
  4%|▍         | frequency blocks : 2/51 [00:00<00:00,   74.27it/s]
  6%|▌         | frequency blocks : 3/51 [00:00<00:00,   68.42it/s]
 10%|▉         | frequency blocks : 5/51 [00:00<00:00,   70.51it/s]
 14%|█▎        | frequency blocks : 7/51 [00:00<00:00,   70.13it/s]
 18%|█▊        | frequency blocks : 9/51 [00:00<00:00,   72.87it/s]
 22%|██▏       | frequency blocks : 11/51 [00:00<00:00,   73.18it/s]
 25%|██▌       | frequency blocks : 13/51 [00:00<00:00,   72.94it/s]
 29%|██▉       | frequency blocks : 15/51 [00:00<00:00,   72.62it/s]
 33%|███▎      | frequency blocks : 17/51 [00:00<00:00,   73.79it/s]
 37%|███▋      | frequency blocks : 19/51 [00:00<00:00,   74.24it/s]
 39%|███▉      | frequency blocks : 20/51 [00:00<00:00,   73.14it/s]
 43%|████▎     | frequency blocks : 22/51 [00:00<00:00,   73.50it/s]
 47%|████▋     | frequency blocks : 24/51 [00:00<00:00,   73.33it/s]
 51%|█████     | frequency blocks : 26/51 [00:00<00:00,   75.29it/s]
 55%|█████▍    | frequency blocks : 28/51 [00:00<00:00,   75.55it/s]
 59%|█████▉    | frequency blocks : 30/51 [00:00<00:00,   75.12it/s]
 63%|██████▎   | frequency blocks : 32/51 [00:00<00:00,   74.99it/s]
 67%|██████▋   | frequency blocks : 34/51 [00:00<00:00,   75.90it/s]
 71%|███████   | frequency blocks : 36/51 [00:00<00:00,   75.22it/s]
 75%|███████▍  | frequency blocks : 38/51 [00:00<00:00,   75.71it/s]
 78%|███████▊  | frequency blocks : 40/51 [00:00<00:00,   75.45it/s]
 82%|████████▏ | frequency blocks : 42/51 [00:00<00:00,   76.32it/s]
 86%|████████▋ | frequency blocks : 44/51 [00:00<00:00,   77.09it/s]
 90%|█████████ | frequency blocks : 46/51 [00:00<00:00,   77.04it/s]
 94%|█████████▍| frequency blocks : 48/51 [00:00<00:00,   76.76it/s]
 98%|█████████▊| frequency blocks : 50/51 [00:00<00:00,   76.42it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00,   75.78it/s]
[Connectivity computation done]

Text(0.5, 0.98, 'GC: [A => B]')

Handling high-dimensional data

An important issue to consider when computing multivariate GC is that the data GC is computed on should not be rank deficient (i.e. must have full rank). More specifically, the autocovariance matrix must not be singular or close to singular.

In the case that your data is not full rank and rank is left as None, an automatic rank computation is performed and an appropriate degree of dimensionality reduction will be enforced. The rank of the data is determined by computing the singular values of the data and finding those within a factor of \(1e^{-6}\) relative to the largest singular value.

Whilst unlikely, there may be scenarios in which this threshold is too lenient. In these cases, you should inspect the singular values of your data to identify an appropriate degree of dimensionality reduction to perform, which you can then specify manually using the rank argument. The code below shows one possible approach for finding an appropriate rank of close-to-singular data with a more conservative threshold.

# gets the singular values of the data
s = np.linalg.svd(raw.get_data(), compute_uv=False)
# finds how many singular values are 'close' to the largest singular value
rank = np.count_nonzero(s >= s[0] * 1e-4)  # 1e-4 is the 'closeness' criteria

Nonethless, even in situations where you specify an appropriate rank, it is not guaranteed that the subsequently-computed autocovariance sequence will retain this non-singularity (this can depend on, e.g. the number of lags). Hence, you may also encounter situations where you have to specify a rank less than that of your data to ensure that the autocovariance sequence is non-singular.

In the above examples, notice how a rank of 5 was given, despite there being 20 channels in the seeds and targets. Attempting to compute GC on the original data would not succeed, given that the resulting autocovariance sequence is singular, as the example below shows.

spectral_connectivity_epochs(
    epochs,
    method=["gc"],
    indices=indices_ab,
    fmin=5,
    fmax=30,
    rank=None,
    gc_n_lags=20,
    verbose=False,
)  # A => B

Traceback (most recent call last):
  File "/home/circleci/project/examples/granger_causality.py", line 415, in <module>
    spectral_connectivity_epochs(
  File "<decorator-gen-312>", line 10, in spectral_connectivity_epochs
  File "/home/circleci/project/mne_connectivity/spectral/epochs.py", line 1383, in spectral_connectivity_epochs
    conn_method.compute_con(indices_use, rank, n_epochs)
  File "/home/circleci/project/mne_connectivity/spectral/epochs_multivariate.py", line 859, in compute_con
    A_f, V = self._autocov_to_full_var(autocov)
  File "/home/circleci/project/mne_connectivity/spectral/epochs_multivariate.py", line 954, in _autocov_to_full_var
    raise RuntimeError(
RuntimeError: the autocovariance matrix is singular; check if your data is rank deficient and specify an appropriate rank argument <= the rank of the seeds and targets

Rigorous checks are implemented to identify any such instances which would otherwise cause the GC computation to produce erroneous results. You can therefore be confident as an end-user that these cases will be caught.

Finally, when comparing GC scores across recordings, it is highly recommended to estimate connectivity from the same number of channels (or equally from the same degree of rank subspace projection) to avoid biases in connectivity estimates. Bias can be avoided by specifying a consistent rank subspace to project to using the rank argument, standardising your connectivity estimates regardless of changes in e.g. the number of channels across recordings. Note that this does not refer to the number of seeds and targets within a connection being identical, rather to the number of seeds and targets across connections.

References

[1] (1,2)

Lionel Barnett and Anil K. Seth. Granger causality for state-space models. Physical Review E, 91(4):040101, 2015. doi:10.1103/PhysRevE.91.040101.

[2] (1,2)

Irene Winkler, Danny Panknin, Daniel Bartz, Klaus-Robert Müller, and Stefan Haufe. Validity of time reversal for testing granger causality. IEEE Transactions on Signal Processing, 64(11):2746–2760, 2016. doi:10.1109/TSP.2016.2531628.

[3]

Peter Whittle. On the fitting of multivariate autoregressions, and the approximate canonical factorization of a spectral density matrix. Biometrika, 50(1-2):129–134, 1963. doi:10.1093/biomet/50.1-2.129.

[4] (1,2)

Stefan Haufe, Vadim V Nikulin, Klaus-Robert Müller, and Guido Nolte. A critical assessment of connectivity measures for eeg data: a simulation study. NeuroImage, 64:120–133, 2013. doi:10.1016/j.neuroimage.2012.09.036.

[5]

Stefan Haufe, Vadim V Nikulin, and Guido Nolte. Alleviating the influence of weak data asymmetries on granger-causal analyses. In Latent Variable Analysis and Signal Separation: 10th International Conference, LVA/ICA 2012, Tel Aviv, Israel, March 12-15, 2012. Proceedings 10, 25–33. Springer, 2012. doi:10.1007/978-3-642-28551-6_4.

[6]

Martin Vinck, Lisanne Huurdeman, Conrado A Bosman, Pascal Fries, Francesco P Battaglia, Cyriel MA Pennartz, and Paul H Tiesinga. How to detect the granger-causal flow direction in the presence of additive noise? NeuroImage, 108:301–318, 2015. doi:10.1016/j.neuroimage.2014.12.017.