Inferring the temporal structure of directed functional connectivity in neural systems: some extensions to Granger causality

TL;DR

This paper extends Granger causality methods to better analyze the temporal dynamics of neural connectivity, enabling more detailed insights into brain information flow across different time scales.

Contribution

It introduces new variants of Granger causality based on multi-step, infinite-future, and single-lag prediction for enhanced temporal analysis.

Findings

New Granger causality variants for temporal analysis

Improved understanding of neural information flow

Systematic approach to neural connectivity timing

Abstract

Neural processes in the brain operate at a range of temporal scales. Granger causality, the most widely-used neuroscientific tool for inference of directed functional connectivity from neurophsyiological data, is traditionally deployed in the form of one-step-ahead prediction regardless of the data sampling rate, and as such yields only limited insight into the temporal structure of the underlying neural processes. We introduce Granger causality variants based on multi-step, infinite-future and single-lag prediction, which facilitate a more detailed and systematic temporal analysis of information flow in the brain.

Tables

Table 1. Table 1: Simple 5 5 5 -variable AR model parameters (Section  5 ).

target ($x$)	source ($y$)	AR lag	AR coefficient
$1$	$2$	$11$	$0.221$
$2$	$1$	$5$	$0.306$
$3$	$1$	$8$	$-0.403$
$4$	$3$	$20$	$-0.215$
$3$	$5$	$4$	$0.352$

Equations

$$\mathbf{P}\left({\mathbf{x}}_{t+1}|{\mathbf{u}}_{-\infty :t}\right)=\mathbf{P}\left({\mathbf{x}}_{t+1}|{\mathbf{u}}_{-\infty :t}^{[y]}\right),$$

$${\mathbf{u}}_t=\sum _{k=0}^{\infty }{B}_k{\mathbf{\epsilon }}_{t-k},\text{ or }{\mathbf{u}}_t=B(z){\mathbf{\epsilon }}_t,$$

$${\mathbf{u}}_t=\sum _{k=1}^{\infty }{A}_k{\mathbf{u}}_{t-k}+{\mathbf{\epsilon }}_t,\text{ or }A(z){\mathbf{u}}_t={\mathbf{\epsilon }}_t,$$

$$\mathbf{E}\left[{\mathbf{u}}_{t+1}|{\mathbf{u}}_{-\infty :t}\right]=\sum _{k=1}^{\infty }{A}_k{\mathbf{u}}_{t+1-k},$$

$${\mathbf{u}}_t^{[y]}=\sum _{k=1}^{\infty }{A}_k^{[y]}{\mathbf{u}}_{t-k}^{[y]}+{\mathbf{\epsilon }}_t^{[y]},\text{ or }{A}^{[y]}(z){\mathbf{u}}_t^{[y]}={\mathbf{\epsilon }}_t^{[y]},$$

$$F_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}=\log \frac{\left|{\Sigma }_{xx}^{[y]}\right|}{\left|{\Sigma }_{xx}\right|}$$

$$F_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}=0⟺A_{xy}(z)\equiv 0.$$

$$H_0:A_{k,xy}=0,k=1,\dots ,p.$$

$$\mathbf{E}\left[{\mathbf{u}}_{t+h}|{\mathbf{u}}_{-\infty :t}\right]=\sum _{k=h}^{\infty }{B}_k{\mathbf{\epsilon }}_{t+h-k},$$

$${\mathbf{\epsilon }}_t^{(h)}=\sum _{k=0}^{h-1}{B}_k{\mathbf{\epsilon }}_{t+h-k}$$

$${\Sigma }^{(h)}=\mathbf{E}\left[{\mathbf{\epsilon }}_t^{(h)}{\mathbf{\epsilon }}_t^{(h)\text{scalebox}0.6\mathrm{T}}\right]=\sum _{k=0}^{h-1}{B}_k\Sigma B_k^{\text{scalebox}0.6\mathrm{T}}.$$

$${\mathbf{u}}_t=\sum _{k=h}^{\infty }{A}_k^{(h)}{\mathbf{u}}_{t-k}+{\mathbf{\epsilon }}_t^{(h)},\text{ or }{A}^{(h)}(z){\mathbf{u}}_t={\mathbf{\epsilon }}_t^{(h)},$$

$$\mathbf{E}\left[{\mathbf{u}}_{t+h}|{\mathbf{u}}_{-\infty :t}\right]=\sum _{k=h}^{\infty }{A}_k^{(h)}{\mathbf{u}}_{t+h-k}.$$

$$A_{h+k}^{(h+1)}=A_{h+k}^{(h)}+A_h^{(h)}{A}_k,h,k=1,2,\dots ,$$

$$F_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}^{(h)}=\log \frac{\left|{\Sigma }_{xx}^{[y](h)}\right|}{\left|{\Sigma }_{xx}^{(h)}\right|},$$

$$F_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}^{(h)}=0⟺A_{xy}^{(h)}(z)\equiv 0.$$

$$F_{\mathbf{y}\to \mathbf{x}}=0⟺F_{\mathbf{y}\to \mathbf{x}}^{(h)}\forall h>0;$$

$$\underset{h\to \infty }{\lim }{F}_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}^{(h)}=0.$$

$$H_0:A_{k,xy}^{(h)}=0,k=h,\dots ,p,$$

$$B_k=A_k+\sum _{\ell =1}^{k-1}{B}_{\ell }{A}_{k-\ell },k=2,3,\dots ,$$

$${\Sigma }_{xx}^{\{h\}}={\left[{\Sigma }^{p,q}\right]}_{xx},p,q=0,\dots ,h-1$$

$${\Sigma }^{p,q}=\sum _{k=0}^{h-1}\sum _{\ell =0}^{h-1}{\delta }_{p-k,q-\ell }{B}_k\Sigma B_{\ell }^{\text{scalebox}0.6\mathrm{T}}.$$

$${\Sigma }_{xx}^{\{h\}}=B_x^{\{h\}}{\Sigma }^{\otimes h}{B}_x^{\{h\}\text{scalebox}0.6\mathrm{T}},$$

$$\begin{gathered} B_x^{\{h\}}= \\ \left[\begin{array}{ccccc}{B}_{0,xu}& 0& 0& \cdots & 0\\ B_{1,xu}& B_{0,xu}& 0& \cdots & 0\\ B_{2,xu}& B_{1,xu}& B_{0,xu}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ B_{h-1,xu}& B_{h-2,xu}& B_{h-3,xu}& \cdots & B_{0,xu}\end{array}\right] \end{gathered}$$

$${\Sigma }^{\otimes h}=\left[\begin{array}{cccc}\Sigma & 0& \cdots & 0\\ 0& \Sigma & \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \Sigma \end{array}\right]$$

$$F_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}^{\{\infty \}}=\underset{h\to \infty }{\lim }{F}_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}^{\{h\}},$$

$$F_{\mathbf{y}\to \mathbf{x}|\mathbf{z}}^{\{h\}}=\log \frac{\left|{\Sigma }_{xx}^{[y]\{h\}}\right|}{\left|{\Sigma }_{xx}^{\{h\}}\right|}.$$

$$H_0:A_{\tau ,xy}=0,$$

$${\mathbf{x}}_t$$

$$+A_{1,xx}{\mathbf{y}}_{t-1}+A_{2,xx}{\mathbf{y}}_{t-2}+\dots +\boxed{A_{\tau ,xy}{\mathbf{y}}_{t-\tau }}+\dots$$

$$+A_{1,xz}{\mathbf{z}}_{t-1}+A_{2,xz}{\mathbf{z}}_{t-2}+{\mathbf{\epsilon }}_{xt}.$$

Full text

1 Introduction

Granger causality (henceforth GC) [21, 22] is a statistical, predictive notion of causal influence originally developed in econometrics, which may be inferred from time-series data, and intuitively interpreted as information flow [1, 2]. Over the past couple of decades it has rapidly become a popular tool for the inference from neurophysiological time series data of time-directed functional (i.e., statistical) relationships in the underlying neural dynamics.

Brain recording modes such as M/EEG, ECoG and fMRI may be characterised as the discrete, regular sampling of continuous-time analogue signals associated with underlying neural processes [6]. Due to variation in biophysical parameters such as axonal length, diameter, conduction velocity, myelination and synaptic delay [31, 9, 11], such processes typically feature signal propagation delays at a range of time scales. Typical application of GC, however, involves prediction only at the time scale of a single time step into the future with respect to the chosen sampling rate. Econometricians have long known that time aggregation of signals engendered by discrete subsampling can induce spurious GC inference [13, 33, 8, 30, 38, 39] and confound detection of actual GCs [3, 36, 46, 6].

There has, in addition, been an awareness that restriction to single-step prediction may obscure the temporal details of causal interactions within a system [25, 27, 16], potentially leading to misinterpretation of GC inferences. Dufour and Renault in particular [16] present a thorough analysis of GC based on multiple prediction horizons, deriving algebraic conditions for (non-)causality between constituent sub-processes at all time scales. Faes et al. [18, 17] present a distinct approach, where causal time scales are explored through causal filtering followed by downsampling. Here, building on [16], we present quantitative statistics for multi-step GC, and demonstrate how they may be estimated in sample and deployed for statistical inference. We also present an infinite-future GC statistic which summarises the total directed connectivity at all accessible predictive time scales, and a single-lag GC which forensically identifies individual causal feedback at specific time lags. The intention is that these statistical tools facilitate unpicking the rich multi-scale temporal details of directed functional interactions in complex neural dynamics.

2 Wiener-Granger causality

Wiener-Granger causality [42, 21, 22] is premised on a notion of causation whereby cause (i) precedes effect, and (ii) contains unique information about the effect. Formally, we suppose given a discrete-time, $n$-dimensional vector111All vector quantities are taken to be column vectors. stochastic process $\bm{u}=\{\left.{\bm{u}_{t}}\,\right|\,{t\in\mathbb{Z}}\}$ representing the “universe of available information”. We introduce the notation $\bm{u}_{t_{1}:t_{2}}$ for the range $\{\left.{\bm{u}_{t}}\,\right|\,{t_{1}\leq t\leq t_{2}}\}$, so that in particular $\bm{u}_{-\infty:t}=\{\left.{\bm{u}_{s}}\,\right|\,{s\leq t}\}$ denotes the infinite history of $\bm{u}$ up to and including time $t$.

Suppose now that $\bm{u}$ is partitioned into non-overlapping sub-processes, $\bm{u}_{t}=[\bm{x}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}\>\bm{y}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}\>\bm{z}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}]^{\mathrm{\scalebox{0.6}{{T}}}}$, of dimension $n_{x},n_{y},n_{z}$ respectively. We say that $\bm{y}$ does not Granger-cause $\bm{x}$ (at time $t$) iff

[TABLE]

where ${\mathbf{P}}\!\left({{\cdot}\,\left|\,{\cdot}\right.}\right)$ denotes conditional distribution, and $\bm{u}^{[y]}=[\bm{x}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}\>\bm{z}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}]^{\mathrm{\scalebox{0.6}{{T}}}}$ the “reduced” universe of information with $\bm{y}$ omitted. Intuitively, (1) says that removing the influence of $\bm{y}$ from the historical information set makes no difference to the statistical distribution of $\bm{x}$ at the next time step, and we say that $\bm{y}$ Granger-causes $\bm{x}$ iff (1) does not obtain. Granger [21, 22] operationalised this definition in terms of (linear) prediction:

$\bm{y}$ Granger-causes $\bm{x}$ iff the history of $\bm{y}$ improves prediction of the future of $\bm{x}$ beyond the extent to which $\bm{x}$ is already predicted by all other available historical information, including that of $\bm{x}$ itself.

Of course in practice, the “universe of available information” will be restricted to a specified set of accessible observables.

Granger causality was subsequently quantified by Geweke [19, 20] as a log-likelihood ratio statistic, and more recently [1, 2] granted an information-theoretic (non-parametric) interpretation in terms of the closely-related transfer entropy [35, 32] (in fact if all stochastic variables are jointly Gaussian [1], the concepts coincide). Under this interpretation—which we prefer—Granger causality represents a “flow of information” from the process $\bm{y}$ to the process $\bm{x}$.

3 Granger-Geweke causality

Suppose now that the process $\bm{u}$ is covariance-stationary, and without loss of generality we assume it to be zero-mean. Then by Wold’s theorem [34, 28], $\bm{u}$ has a moving average (MA) representation

[TABLE]

where $\bm{\varepsilon}$ is a white noise process with nonsingular covariance matrix $\Sigma={\mathbf{E}}\!\left[{\bm{\varepsilon}_{t}\bm{\varepsilon}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}}\right]$, $z$ is the lag (backshift) operator222Note that in the literature, the lag operator is sometimes taken as $z^{-1}$. In the spectral domain, $z$ may be viewed as residing on the unit circle in the complex plane: $z=e^{-i\omega}$, where $\omega$ is the phase angle in radians. (so that $z\cdot\bm{\varepsilon}_{t}=\bm{\varepsilon}_{t-1}$, etc.), and the MA operator (transfer function) is given by $B(z)=\sum_{k=0}^{\infty}B_{k}z^{k}$ with $B_{k}$ the MA coefficient matrices, and $B_{0}=I$ (the identity matrix), so that $B(z)$ is causal (does not reference the future). We also assume the minimum-phase condition that $B(z)$ is nonsingular on the closed unit disc in the complex plane, so that the MA representation (2) may be inverted to yield a stable, causal autoregressive (AR) representation

[TABLE]

where $A(z)=B(z)^{-1}=I-\sum_{k=0}^{\infty}A_{k}z^{k}$ is also nonsingular on the closed unit disc (see, e.g., [29, 34, 19]).

Granger considered prediction in the linear least-squares sense. The optimal linear prediction of $\bm{u}_{t+1}$ given its history $\bm{u}_{-\infty:t}$ is the conditional expectation [23]

[TABLE]

with residual prediction error $\bm{\varepsilon}_{t+1}$. Following [19], prediction error is quantified by the determinant $\left|{\Sigma}\right|$ of the residuals covariance matrix, also known as the generalised variance [43, 7]. Considering now the partition $\bm{u}_{t}=[\bm{x}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}\>\bm{y}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}\;\bm{z}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}]^{\mathrm{\scalebox{0.6}{{T}}}}$, the optimal linear prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t+1}}\,\left|\,{\bm{u}_{-\infty:t}}\right.}\right]$ of $\bm{x}_{t+1}$ given the full history $\bm{u}_{-\infty:t}$ has prediction error $\bm{\varepsilon}_{x,t+1}$ with generalised variance $\left|{\Sigma_{xx}}\right|$ (here subscript ‘$x$’ denotes $x$-component). We contrast this with the optimal prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t+1}}\,\left|\,{\bm{u}^{[y]}_{-\infty:t}}\right.}\right]$ of $\bm{x}_{t+1}$ on the reduced universe of historical information $\bm{u}^{[y]}=[\bm{x}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}\>\bm{z}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}]^{\mathrm{\scalebox{0.6}{{T}}}}$, where $\bm{y}$ is omitted; cf., (1) . This is derived from the reduced AR representation

[TABLE]

so that the optimal prediction of $\bm{u}^{[y]}_{t+1}$ on its own history is ${\mathbf{E}}\!\left[{{\bm{u}^{[y]}_{t+1}}\,\left|\,{\bm{u}^{[y]}_{-\infty:t}}\right.}\right]=\sum_{k=1}^{\infty}A^{[y]}_{k}\bm{u}^{[y]}_{t+1-k}$, and the generalised variance for the optimal prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t+1}}\,\left|\,{\bm{u}^{[y]}_{-\infty:t}}\right.}\right]$ is $\left|{\Sigma^{[y]}_{xx}}\right|$ where $\Sigma^{[y]}={\mathbf{E}}\!\left[{\bm{\varepsilon}^{[y]}_{t}\bm{\varepsilon}^{{[y]}\mathrm{\scalebox{0.6}{{T}}}}_{t}}\right]$. Following [20] the Granger-Geweke causality statistic is defined as

[TABLE]

and we have, in particular,

[TABLE]

In finite sample, where the infinite histories are truncated at some model order $p$ and models (3) and (5) estimated by maximum likelihood (e.g., an OLS), the estimated generalised variances are proportional to the likelihoods, and the sample estimator $\hat{F}_{\bm{y}\to\bm{x}|\bm{z}}$ is a log-likelihood ratio statistic. In this scenario, (3, 5) are nested linear autoregression models, and the null hypothesis of vanishing Granger causality is

[TABLE]

Thus, by the standard large-sample theory [44, 40], under the null hypothesis (8) the maximum-likelihood estimator $\hat{F}_{\bm{y}\to\bm{x}|\bm{z}}$, scaled by sample size, converges in distribution to a central $\chi^{2}(d)$ with degrees of freedom given by $d=pn_{x}n_{y}$, and a non-central $\chi^{2}(d;\lambda)$ with non-centrality parameter $\lambda=F_{\bm{y}\to\bm{x}|\bm{z}}$ under the alternative hypothesis333We remark that in sample, the scaled lack-of-fit sum of squares $[\operatorname{trace}\!\left({\Sigma^{[y]}_{xx}}\right)-\operatorname{trace}\!\left({\Sigma_{xx}}\right)]/\operatorname{trace}\!\left({\Sigma_{xx}}\right)$ is asymptotically $F$-distributed under $H_{0}$, furnishing an alternative and more statistically powerful test for the null. However, the F-statistic lacks an information-theoretic interpretation, as well as some crucial invariance properties [7, 3] of the log-likelihood ratio form, and is thus less satisfactory as a measure of magnitude of Granger-causal effect.. We note that $F_{\bm{y}\to\bm{x}|\bm{z}}$ is strictly non-negative, and thus biased in sample.

The above finite-sample analysis assumes that the AR models (3, 5) are independently estimated. It is, however, known [12] that this may be problematic, in particular for spectral (frequency-domain) Granger-Geweke causality [19, 20] (which we do not address here). In fact, from the Spectral Factorisation Theorem [29] it follows that the reduced model (5) may be deduced from the full model (3), leading to more powerful and less biased GC estimators. There are several approaches to effecting this computationally, in the frequency domain [45, 14, 15] and in the time domain [41, 4]. More recently, [5, 39] show how this may be efficiently accomplished using state-space methods [24]444Sample statistics derived from single full-model estimation with spectral factorisation, however, fail to satisfy the requirements for the large-sample theory; in lieu of known distributions for these estimators, independent estimates of the full and reduced models or standard subsampling/surrogate methods, may be considered preferable for statistical inference.; there are, furthermore, other compelling reasons to estimate Granger-Geweke causality via state-space rather than AR modelling [5, 39].

4 Multi-step Granger causality

In the traditional approach, Granger causality is usually considered only in terms of one-step-ahead prediction [cf., (1, 4)]; but see, e.g., [27, 16]. However, as noted in Section 1 the duration of a single time step will vary according to the sampling rate, and the magnitude of reported Granger-Geweke causality will depend crucially on the relationship between sampling frequency and underlying time scales of neural signal transmission [6]. This suggests we examine more closely Granger causality based on an arbitrary future prediction horizon. A notion of (non-)causality consonant with the measure we consider was introduced in [16]; here, for the first time (as far as we are aware), we quantify this notion with a Granger-Geweke statistic.

We require an expression for ${\mathbf{E}}\!\left[{{\bm{u}_{t+h}}\,\left|\,{\bm{u}_{-\infty:t}}\right.}\right]$, $h=1,2,\ldots$; that is, optimal linear prediction at an arbitrary future prediction horizon $h$ (but note that the historical predictor set $\bm{u}_{-\infty:t}$ remains the same as for conventional $1$-step GC). In this case the $h$-step optimal prediction is more simply expressed in terms of the MA, rather than AR representation [cf., (4)]. We have [23]

[TABLE]

with residual errors

[TABLE]

[henceforth we use the round-bracket ‘${\bm{(}h\bm{)}}$’ to indicate a prediction horizon $h$ steps into the future]. Note that in general $\bm{\varepsilon}^{\bm{(}h\bm{)}}$ will not be a white noise process. The residuals covariance matrix is given by

[TABLE]

Setting $B^{\bm{(}h\bm{)}}(z)=\sum_{k=0}^{h-1}B_{k}z^{k}$ and $A^{\bm{(}h\bm{)}}(z)=B^{\bm{(}h\bm{)}}(z)A(z)=I-\sum_{k=h}^{\infty}A^{\bm{(}h\bm{)}}_{k}z^{k}$, we may derive the $h$-lagged AR form [cf., (3)]

[TABLE]

and the AR expression for the optimal $h$-step linear prediction [cf., (4)]

[TABLE]

The $A^{\bm{(}h\bm{)}}_{k}$ satisfy the recursion relations [16]

[TABLE]

with $A^{\bm{(}1\bm{)}}_{k}=A_{k}$.

We now define $h$-step Granger-Geweke causality by analogy with (6) as [6]

[TABLE]

where $\Sigma^{{[y]}{\bm{(}h\bm{)}}}={\mathbf{E}}\!\left[{\bm{\varepsilon}^{{[y]}{\bm{(}h\bm{)}}}_{t}\bm{\varepsilon}^{{[y]}{\bm{(}h\bm{)}}\mathrm{\scalebox{0.6}{{T}}}}_{t}}\right]=\sum_{k=0}^{h-1}B^{[y]}_{k}\Sigma^{[y]}B^{{[y]}\mathrm{\scalebox{0.6}{{T}}}}_{k}$, and we have [cf., (7)]

[TABLE]

In contrast to the $1$-step case (7), this condition will generally be nonlinear—specifically, a series of matrix polynomial identities of order $h$—in the AR coefficients $A_{k}$. In the unconditional case $\bm{z}=\emptyset$, it may be shown [37, 10] that

[TABLE]

however, in the conditional case, neither implication holds in general [27, 16]. We note also from (2, 10, 11), that as $h\to\infty$, both $\Sigma^{\bm{(}h\bm{)}}_{xx}$ and $\Sigma^{{[y]}{\bm{(}h\bm{)}}}_{xx}\to{\mathbf{E}}\!\left[{\bm{x}_{t}\bm{x}_{t}^{\mathrm{\scalebox{0.6}{{T}}}}}\right]$, the covariance matrix of $\bm{x}$ itself, implying [6]

[TABLE]

Related analysis in a continuous-time scenario [6] suggests that convergence in (18) is exponential.

From (12) it follows that we again have nested ($h$-step AR) models, so that in finite sample with truncation at $p\geq h$, the null hypothesis of vanishing $h$-step Granger causality is [cf., (8)]

[TABLE]

and again the scaled maximum-likelihood sample estimator for $\hat{F}_{\bm{y}\to\bm{x}}^{\bm{(}h\bm{)}}$ will be asymptotically $\chi^{2}(d)$ under the null hypothesis (19), now with $d=(p-h+1)n_{x}n_{y}$. Computationally, multi-step GC may be estimated from AR or state-space models, using (11, 15). For AR modelling, the MA coefficients may be calculated recursively using

[TABLE]

with $B_{1}=A_{1}$. For state-space modelling, calculation of the $B_{k}$ is even more straightforward (see [5], eq. 4).

5 Full-future Granger causality

Historically, the main emphasis of Granger causality analysis, especially in the econometrics literature, has been on statistical inference of (non-)causality. However, in light of the more recent interpretation of GC as a measure of information flow [1, 2], the Granger-Geweke statistic stands as an effect size, which quantifies this information flow. This perspective seems to us particularly appropriate and intuitive with regard to functional analysis of neural systems. The conventional $1$-step prediction GC statistic, however, may be considered potentially misleading as a comparative effect size, insofar as it fails to take into account neural time scales and their interplay with sampling rate. It would thus be useful to have (in addition to the multi-step GC of Section 4), a summary GC measure of the total information flow between variables; i.e., from infinite past to infinite future. This motivates our introduction of a “full-future” GC measure, based on past-conditional prediction of the infinite future; that is, ${\mathbf{E}}\!\left[{{\bm{u}_{t+1:\infty}}\,\left|\,{\bm{u}_{-\infty:t}}\right.}\right]$.

We may calculate that the residuals covariance matrix of the prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t+1:t+h}}\,\left|\,{\bm{u}_{-\infty:t}}\right.}\right]$ of the future of $\bm{x}$ up to horizon $t+h$ from the full process history $\bm{u}_{-\infty:t}$, is given by the ($h\times h$)-block matrix

[TABLE]

[note: we use curly braces $\bm{\{}h\bm{\}}$ to distinguish the full-future prediction horizon from the multi-step horizon $\bm{(}h\bm{)}$], where

[TABLE]

This may be written

[TABLE]

with

[TABLE]

where the index $u$ denotes all components $\{x,y,z\}$, and

[TABLE]

is ($h\times h$)-block-diagonal.

For the reduced prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t+1:t+h}}\,\left|\,{\bm{u}^{[y]}_{-\infty:t}}\right.}\right]$, we obtain $\Sigma^{{[y]}{\bm{\{}h\bm{\}}}}_{xx}$ the same way, replacing $\Sigma$ with $\Sigma^{[y]}$ and $B_{k}$ with $B^{[y]}_{k}$, and we define the full-future Granger causality as

[TABLE]

where

[TABLE]

We conjecture that under appropriate conditions (cf., Section 2) the limit in (27) exists; this has been verified by extensive simulation.

In sample with finite AR model order, the null hypothesis for vanishing $F^{\bm{\{}\infty\bm{\}}}{\bm{y}\to\bm{x}|\bm{z}}$ is identical to the null (8) for $1$-step Granger causality (6). This follows from the recursion relations (14) and expanding out the prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t+1:t+h}}\,\left|\,{\bm{u}_{-\infty:t}}\right.}\right]$. Thus the statistic is not useful in its own right for statistical inference, and should rather be considered an informative quantitative measure of total past $\,\to\,$ future information flow between two variables.

We have not found a closed formula for the determinants in (27), but they may be approximated numerically; extensive simulations suggest that, although the size of the matrices $B^{\bm{\{}h\bm{\}}}_{x}$ scale quadratically in $h$, convergence to the limit in (27) is again exponential (cf., Fig. 2). We remark that in general, $F^{\bm{\{}\infty\bm{\}}}_{\bm{y}\to\bm{x}|\bm{z}}$ will not be equal to the sum $\sum_{h=1}^{\infty}F^{\bm{(}h\bm{)}}_{\bm{y}\to\bm{x}|\bm{z}}$ of multi-step Granger causalities, since the residuals $\bm{\varepsilon}^{\bm{(}h\bm{)}}_{t}$ of the latter (10) for different $h$ will in general be correlated, so that $\left|{\Sigma^{\bm{\{}h\bm{\}}}}\right|\neq\prod_{k=1}^{h}\left|{\Sigma^{\bm{(}k\bm{)}}}\right|$.

We demonstrate multi-step and full-future GC with a simple AR model with $n=5$ variables (Fig. 1),

and model order $p=20$. All AR coefficients were zero except for the lag-$1$ self-regression terms $A_{1,ii}$, and the lagged coefficients set out in TABLE 1.

The only non-zero multi-step and full-future GCs are plotted in Fig. 2, calculated according to (15) and (27) respectively for prediction horizon $h=1,\ldots,32$, with, for each directed pair of variables $x,y$, full conditioning on all remaining variables $\bm{z}$. Note that the both the multi-step and full-future GCs coincide with the conventional $1$-step GC (6) at prediction horizon $h=1$. Vertical grey lines indicate the AR lag of the causal interaction (TABLE 1). We see that (cf., [6]) the $F^{\bm{(}h\bm{)}}_{y\to x|\bm{z}}$ decay rapidly to zero beyond the causal horizon, while the $F^{\bm{\{}h\bm{\}}}_{y\to x|\bm{z}}$ rise and then quickly plateau beyond the causal horizon to the limiting value $F^{\bm{\{}\infty\bm{\}}}_{y\to x|\bm{z}}$ (26).

6 Single-lag Granger causality

A more fine-grained analysis of directed functional connectivity may be interrogated as follows: if a variable $\bm{y}$ Granger-causes the variable $\bm{x}$, at which specific time scale(s) is causal feedback concentrated? Note that multi-step GC (Section 4) does not directly address this question, since there all lags of a predictive (source) variable are considered together. Rather, for a specific lag $\tau>0$, we consider a Granger statistic based on the null hypothesis

[TABLE]

where $A_{k}$ are the AR coefficient matrices in (3). We thus compare the optimum prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t}}\,\left|\,{\bm{u}_{-\infty:t-1}}\right.}\right]$ with the optimum prediction ${\mathbf{E}}\!\left[{{\bm{x}_{t}}\,\left|\,{\bm{u}^{\bm{<}y;\tau\bm{>}}_{-\infty:t-1}}\right.}\right]$, where the superscript ‘$\bm{<}\!y;\tau\!\bm{>}$’ indicates that the single lag $\bm{y}_{t-\tau}$ of $\bm{y}$ is omitted from the historical predictor set $\bm{u}_{-\infty:t-1}$. To make this clearer, consider the $x$-component of the AR representation (3) of $\bm{u}_{t}$:

[TABLE]

The reduced AR representation then omits the boxed lag-$\tau$ $\bm{y}$ regressor, and we define the single-lag Granger causality as

[TABLE]

where $\Sigma^{\bm{<}y;\tau\bm{>}}_{xx}$ is the residuals covariance matrix for the reduced AR model. We note that $F^{\bm{<}\!\tau\!\bm{>}}_{\bm{y}\to\bm{x}|\bm{z}}=0\;\forall\tau>0\iff F_{\bm{y}\to\bm{x}|\bm{z}}=0$

The regression (29) with the null condition (28) represents a nested linear model, so that the large-sample theory applies, and the scaled sample estimator $\hat{F}^{<\!\tau\!>}_{\bm{y}\to\bm{x}|\bm{z}}$ will thus be asymptotically $\chi^{2}(d)$ with $d=n_{x}n_{y}$. We remark that interpretation of $F^{\bm{<}\!\tau\!\bm{>}}_{\bm{y}\to\bm{x}|\bm{z}}$ as an effect size for a putative “information flow” is somewhat moot; we may prefer to consider $F^{\bm{<}\!\tau\!\bm{>}}_{\bm{y}\to\bm{x}|\bm{z}}$ simply as a test statistic for inference of (the absence of) a causal feedback from source to target variable at the given lag.

Unlike the previous GC measures, we do not have (given full-model parameters) a construction for a state-space model which represents the reduced model (28). The reduced model parameters may, however, still be solved computationally from the Yule-Walker equations [27]. The full-model Yule-Walker equations up to lag $q$ yield

[TABLE]

with $\Gamma_{k}={\mathbf{E}}\!\left[{\bm{u}_{t}\bm{u}_{t-k}^{\mathrm{\scalebox{0.6}{{T}}}}}\right]$, $k=\ldots,-2,-1,0,1,2,\ldots$ the autocovariance sequence—which may itself be derived from the (estimated) full-model AR coefficients [4]—and

[TABLE]

The reduced Yule-Walker solution for $\Sigma^{\bm{<}y;\tau\bm{>}}$ is then obtained as per (31), after deleting the $y$-columns of the $\tau$-th block-column in $\bm{\Gamma}_{q}$, and the $y$-rows/columns of the $(\tau-1)$-th block-row/column in $\bm{\Lambda}_{q}$. Even though $\bm{\Lambda}_{q}$ may be quite large555For reasonable numerical precision we need sufficient lags $q$ that $\Gamma_{k}\approx 0$ for $k>q$, which will in turn depend on the spectral radius of $A(z)$ [26]; see e.g., [4]., it is positive-definite Toeplitz and thus may be Cholesky-decomposed and efficiently inverted.

We envisage estimating $F^{\bm{<}\!\tau\!\bm{>}}_{\bm{y}\to\bm{x}|\bm{z}}$ from the data for $\tau=1,\ldots,p$ in turn (where the maximum lag $p$—the model order for the full AR model (29)—is selected via a standard scheme), in order to ascertain the time scale(s) at which $\bm{y}$ influences $\bm{x}$. See Fig. 3, where the $F^{\bm{<}\!\tau\!\bm{>}}_{y\to x|\bm{z}}$, $x,y=1,\ldots,n$, $x\neq y$, are estimated in sample for a data sequence of length $1000$ generated from the AR model of Section 5.

Here $\bm{z}$ denotes all other variables except the given $x,y$, so that every directed pairwise GC is conditioned on all remaining variables, yielding the “Granger-causal graph” [4] at all lags up to $p=20$. Likelihood-ratio single-lag GC statistics (blue boxes) were calculated for separate OLS estimates of the full and (for each $j,\tau$) reduced models (29) using the (known) model order $p=20$, while analytic GCs for the model (black horizontal bars) were calculated from the actual model parameters (TABLE 1) using the Yule-Walker procedure described above with $q=175$ autocovariance lags, which was sufficient to ensure that the $\Gamma_{k}$ decay to near-machine precision. The red horizontal lines mark the critical GC level for rejection of the null hypotheses (28) of zero single-lag GC at significance $\alpha=0.05$ according to the $\chi^{2}(1)$ estimator distribution, assuming a Bonferroni correction for all $pn(n-1)$ hypotheses. We see that statistical inference of the $F^{<\!\tau\!>}_{y\to x|\bm{z}}$ correctly identifies the causal lags as well as directed functional connectivity in the model (Fig. 1).

7 Conclusions

In this article we address the issue of how, in an empirical scenario, we may go beyond Granger-causal inference restricted to the time scale prescribed by the data sampling rate, to obtain a more detailed picture of Granger-causal interactions at multiple times scales underlying the measured neurophysiological process. Thus our multi-step measure $F^{\bm{(}h\bm{)}}_{\bm{y}\to\bm{x}|\bm{z}}$ reflects information flow between variables at a specific future time horizon $h$, while the single-lag measure $F^{\bm{<}\!\tau\!\bm{>}}_{\bm{y}\to\bm{x}|\bm{z}}$ identifies the precise time lag(s) at which a specific Granger-causal interaction operates. Via a simple didactic model, we demonstrate, respectively, how the underlying time scales are reflected via the multi-step statistic (Fig. 2), and may be explicitly inferred from the data (Fig. 3). In addition, in our full-future GC measure $F^{\bm{\{}\infty\bm{\}}}_{\bm{y}\to\bm{x}|\bm{z}}$, we present a useful summary measure of effect size for the total past $\to$ future information flow between variables. All measures are fully conditioned on (accessible) exogenous variables, so that only direct functional relationships are reported. We describe how our measures may be estimated computationally from time-series data, and (where appropriate) their asymptotic sampling distributions. We propose these measures as useful additions to the directed functional analysis toolbox, insofar as they stand to elucidate time-dependant causal interactions in neurophysiological processes of interest to neuroscientific research. We encourage future research into the behaviour of our measures for more realistic data, where we should expect causal interactions at multiple distributed lags [6].

Acknowledgments

We are grateful to the Dr. Mortimer and Theresa Sackler Foundation, which supports the Sackler Centre for Consciousness Science. The authors would also like to thank Stefan Haufe for useful discussions.