Solved problems and remaining challenges for Granger causality analysis in neuroscience: A response to Stokes and Purdon (2017)
Lionel Barnett Adam B. Barrett Anil K. Seth


































Sackler Centre for Consciousness Science and Department of Informatics, 

University of Sussex, Brighton BN1 9QJ, UK
a.k.seth@sussex.ac.uk (correspondence)


Granger-Geweke causality (GGC) is a powerful and popular method for identifying directed functional (`causal') connectivity in neuroscience. In a recent paper, <cit.> raise several concerns about its use. They make two primary claims: (1) that GGC estimates may be severely biased or of high variance, and (2) that GGC fails to reveal the full structural/causal mechanisms of a system. However, these claims rest, respectively, on an incomplete evaluation of the literature, and a misconception about what GGC can be said to measure. Here we explain how existing approaches [as implemented, for example, in our popular MVGC software <cit.>] resolve the first issue, and discuss the frequently-misunderstood distinction between functional and effective neural connectivity which underlies Stokes and Purdon's second claim.



Granger causality functional connectivity effective connectivity statistical inference













Granger-Geweke causality (GGC) is a powerful analysis method for inferring directed functional (`causal') connectivity from time-series data, which has become increasingly popular in a variety of neuroimaging contexts <cit.>.  GGC operationalises a statistical, predictive notion of causality in which causes precede, and help predict their effects. When implemented using autoregressive modelling, GGC can be computed in both time and frequency domains, in both bivariate and multivariate (conditional) formulations. Despite its popularity and power, the use of GGC in neuroscience and neuroimaging has remained controversial. In a recent paper, <cit.> raise two primary concerns: (1) that GGC estimates may be severely biased or of high variance, and (2) that GGC fails to reveal the full structural/causal mechanisms of a system. Here, we explain why these concerns are misplaced.

Regarding the first claim, <cit.> describe how bias and variance in GGC estimation arise from the use of separate, independent, full and reduced regressions. However, this problem has long been recognised <cit.> and, moreover, has already been solved by methods which derive GGC from a single full regression[We note here that the “partition matrix” solution proposed by <cit.> is incorrect; see, , <cit.>.]. These methods essentially extract reduced model parameters from the full model via factorisation of the spectral density matrix. Well-documented approaches include Wilson's frequency-domain algorithm <cit.>, Whittle's time-domain algorithm <cit.>, and a state-space approach which devolves to solution of a discrete-time algebraic Riccati equation <cit.>. Thus, the source of bias and variance discussed in <cit.> has already been resolved.

This is clearly illustrated in fig:sgc, where we plot estimated frequency-domain GGC for the 3-node VAR model in <cit.>, Example 1, using the single-regression state-space method <cit.>. We remark that identical results are obtained using the time-domain spectral factorisation method of <cit.>, as implemented in the current (v1.0, 2012) release of the associated MVGC Matlab^ software package <cit.>. fig:sgc may be directly compared with Fig. 2 in <cit.>; we see clearly that all estimates are strictly non-negative, and that exaggerated bias and variance associated with the dual-regression approach are absent. Therefore,  <cit.> are in error when they state that “Barnett and Seth […] have proposed fitting the reduced model and using it to directly compute the spectral components …”. This is important to note because our MVGC toolbox has been widely adopted within the community, with > 3,500 downloads and a significant number of high-impact research publications using the method <cit.>. Thus, we can reassure users of the toolbox that problems of bias and variance as described by <cit.> do not apply.


Sample variance is, of course, still evident, as is bias due to non-negativity of the GGC sample statistic (which may be countered by standard surrogate data methods), but both remain well below their minimum values across all model orders for the dual-regression case <cit.>. fig:bivar further compares bias and variance of time-domain GGC for the example system for single and dual regressions, at model order 3, across a wide-range of time-series lengths. A single regression consistently leads to substantially less bias and variance, except at high time-series lengths, where there is a drop off of bias and variance for both methods.



<cit.> do correctly identify a fundamental cause of the problem with dual-regression GGC estimation: even if the full process is a finite-order autoregression, the reduced process will generally not be finite-order autoregressive; rather, it will be VARMA, or equivalently, a finite-order state-space process <cit.> – which may be poorly modelled as a finite-order VAR <cit.>. The problem is in fact more pervasive than this: the full process itself may have a strong moving-average (MA) component and be poorly-modelled as a finite-order VAR. This is because common features of neurophysiological data acquisition, sampling and preprocessing procedures such as subsampling and other temporal aggregation, filtering, measurement noise and sub-process extraction will all, in general, induce an MA component <cit.>. This is particularly pertinent to fMRI data, where the haemodynamic response acts as a slow, MA filter. Fortunately, the state-space and non-parametric approaches mentioned above handle VARMA data parsimoniously, hence avoiding this problem.

Moving on to the second claim, Stokes and Purdon note that GGC reflects a combination of `transmitter' and `channel' dynamics, and is independent of `receiver' dynamics. Again, this independence has been previously identified, as a direct consequence of the invariance of GGC under certain affine transformations <cit.>. But why should this independence matter? They suggest that it runs “counter to intuitive notions of causality intended to explain observed effects” since, according to them, “neuroscientists seek to determine the mechanisms that produce `effects' within a neural system or circuit as a function of inputs or `causes' observed at other locations”.  In fact, this view resonates more strongly with approaches such as Dynamic Causal Modelling <cit.>—usually characterised as `effective connectivity'—which attempt to find the optimal mechanistic (circuit-level) description that explains observed data. GGC, on the other hand, models dependencies among observed responses and is therefore an example of (directed) `functional connectivity' <cit.>. Essentially, the distinction is between making inferences about an underlying physical causal mechanism <cit.> and making inferences about directed information flow <cit.>. DCM is able to deliver evidence for circuit-level descriptions of neural mechanism from a limited repertoire of tightly-framed hypotheses, which must be independently motivated and validated <cit.>; it is, in particular, unsuited to exploratory analyses. GGC, on the other hand, is data-driven and “data-agnostic” (it makes few assumptions about the generative process, beyond that it be reasonably parsimoniously modelled as a linear stochastic system), and as such is well-suited to exploratory analyses. It delivers an information-theoretic interpretation of the neural process which is both amenable to statistical inference, and which also stands as an effect size for directed information flow between components of the system <cit.>. Both approaches address valid questions of interest to neuroscientific analysis.

Concluding, GGC represents a conceptually satisfying and statistically powerful method for (directed) functional connectivity analysis in neuroscience and neuroimaging. Currently available implementations [, <cit.>] deal appropriately with issues of bias and variance associated with dual regression methods. However, a range of additional challenges remain in further developing this useful technique. These include issues of stationarity, linearity and exogenous influences, as noted by <cit.>, and in addition the influences of noise, sampling rates and temporal/spatial aggregation engendered by neural data acquisition <cit.>.



§ ACKNOWLEDGEMENTS


ABB is funded by EPSRC grant EP/L005131/1. All authors are grateful to the Dr. Mortimer and Theresa Sackler Foundation, which supports the Sackler Centre for Consciousness Science.


model2-names












































































