Reliability of the Granger causality inference

How to characterize information flows in physical, biological, and social systems remains a major theoretical challenge. Granger causality (GC) analysis has been widely used to investigate information flow through causal interactions. We address one of the central questions in GC analysis, that is, the reliability of the GC evaluation and its implications for the causal structures extracted by this analysis. Our work reveals that the manner in which a continuous dynamical process is projected or coarse-grained to a discrete process has a profound impact on the reliability of the GC inference, and different sampling may potentially yield completely opposite inferences. This inference hazard is present for both linear and nonlinear processes. We emphasize that there is a hazard of reaching incorrect conclusions about network topologies, even including statistical (such as small-world or scale-free) properties of the networks, when GC analysis is blindly applied to infer the network topology. We demonstrate this using a small-world network for which a drastic loss of small-world attributes occurs in the reconstructed network using the standard GC approach. We further show how to resolve the paradox that the GC analysis seemingly becomes less reliable when more information is incorporated using finer and finer sampling. Finally, we present strategies to overcome these inference artifacts in order to obtain a reliable GC result.


Introduction
Information flow plays a central role in physical and biological systems, such as in gene regulation and cortical computation in the brain. How to characterize information propagation within dynamical systems is a great theoretical challenge. There has been a long history of studying causal relations between time series [1][2][3][4][5][6]. Recently, Granger causality (GC) analysis has proven itself as a powerful method, widely used for extracting causal connectivity from data in fields ranging from physics [7][8][9][10] and biology [11][12][13][14], to social science [15][16][17][18].
In general, it is often quite involved to acquire reliable GC inference. Usually the GC framework requires the time series to be linear. For linear systems, the directed GC corresponds to the structural connectivity that mediates physical interactions within a system. For nonlinear systems, the relationship between the causal connectivity and the structural connectivity is still under active investigation. Extending the GC theory to nonlinear, non-Gaussian times series remains a central theoretical endeavor [19][20][21]. Note that the notion of GC is statistical rather than structural, i.e. it identifies directed statistical causal connectivity through statistical features of time series. As a statistical method, one of the important issues is how to collect time series for dynamical variables in order to faithfully capture the information flow within the dynamics. In this work, we study the reliability of GC inference for both linear and nonlinear systems and demonstrate that there are inference hazards in the GC analysis arising from the manner in which a continuous dynamical process is projected to a discrete process.
Because most dynamical quantities are continuous in time and GC values are generally evaluated using discrete time series sampled from these continuous time processes, we investigate the reliability of GC evaluation by studying the GC value as a function of the sampling interval length τ (we will term this function the GC sampling structure). We show that, for both linear and nonlinear dynamics, there are surprisingly common features in the GC sampling structure: (i) when by the design of our system there is a causal flow, yet the GC value can vanish over a certain set of τ, yielding an incorrect causal inference; or (ii) conversely, when there is no causal influence by construction, the GC value can become of appreciable size for some ranges of τ, yielding again a potentially erroneous causal inference. Clearly, these phenomena greatly complicate the interpretation of the GC inference and potentially produce opposing causality conclusions when using different sampling τʼs for empirical data. These issues have yet to receive sufficient attention despite the wide use of GC analysis in many fields. It is important to examine these issues in a GC analysis to avoid erroneous causal inference. We use a small-world network reconstruction as an example to illustrate this hazard. As is well known, many important biological or social networks possess small-world attributes [22,23]. The small-world topology is generally linked to efficient parallel information processing on both local and global levels. We show that an incorrect inference can arise, i.e. the small-world attributes are completely lost in the reconstructed network when GC analysis is blindly applied to the network topology reconstruction. Clearly, these sampling issues should also arise in data processing for GC evaluations using low-pass filtering or down-sampling for many types of imaging data, such as the BOLD signal in functional magnetic resonance imaging (fMRI) [24,25].
To remove the inference artifacts in the GC analysis, it seems that a natural resolution is to incorporate as much information as possible by using sufficiently fine sampling. However, as we show below, this approach produces a paradox, i.e. a finer sampling may not imply a more reliable GC inference. In what follows, we will use idealized models to study the mechanisms underlying these phenomena and discuss an approach that can eliminate these artifacts to obtain a reliable GC relation.

GC sampling structures
GC aims to analyze the causal influence of one time series, X t , on another, Y t , by examining whether the statistical prediction of Y t can be improved when the past information of X t is incorporated into it [2,26,27]. First, to fix the notation, we briefly recapitulate the GC analysis [2,27,28]. For two stationary zero-mean discrete time series { } X t and { } Y t , a statistical prediction can usually be accomplished using the autoregression where ϵ t 2 and η t 2 are residuals representing the prediction errors after we consider the shared history for both time series. The sequences ϵ 2 , i.e. the prediction of one time series cannot become worse (i.e. have greater residual variance) after information from the other time series is incorporated, The total To illustrate the application of GC in nonlinear dynamics, we choose network dynamics of conductance-based, pulse-coupled integrate-and-fire (I&F) neurons [29] where v i and g i are the membrane potential and conductance for the ith neuron, respectively, g L is the leak conductance, ϵ L and ϵ E are the resting and reversal potentials, respectively, and σ is the conductance time constant. When v i ϵ ∈ ( ) V , L th , a neuron evolves according to equation (2). When v i reaches the threshold V th , it is reset to the resting potential ϵ L for a refractory period τ ref . . We demonstrate below how we can reliably assess this consistency by first addressing artifacts arising from the sampling in the evaluation of GC values.
Note that the discrete time series used in the GC analysis is usually sampled from continuous-time processes. The sampling effects can be captured by the GC sampling structure, i.e. the GC values, → F y x , → F x y , as functions of the sampling interval τ, denoted by τ For the nonlinear network dynamics of two neurons, from its GC sampling structure as shown in figure 1(a), we can observe the following properties: (I) GC is an oscillatory, decaying function of τ, (II) GC becomes nearly zero, or zero, periodically.
From this GC sampling structure, clearly, one is confronted with the question of what a proper τ is to obtain the GC value for causal inference or network topology extraction.
Properties I-II are general, not limited to nonlinear dynamics, and they also appear in many of the linear models we have examined. Figure 1(b) illustrates such an example of linear (second-order regressive) models: where ε t and η t are Gaussian white noise with zero mean, and Here, 'sampling interval length' τ = ( ) k refers to a 'coarse-grained' discrete time series constructed using an original datum point after skipping − k 1 points in between. As shown in figure 1(b), clearly, the model also possesses properties I-II.
Next, we construct idealized models for which we have perfect control of causal influence and use them to ascertain the phenomena observed in the GC sampling structure analytically. Using the lag operator L, i.e.
where ϵ t and η t are independent white-noise series, so that Clearly, the structural connectivity in the system (4) means that Y t causally influences X t and is not influenced by X t . We first consider model I: , which are GC computed using the time series sampled with interval length k. Therefore, for the sampling interval k, we have for large k, which clearly possesses properties I-II. In particular, for suitable values of β and ϕ, the GC value vanishes on a certain set of k. Figure 1(c) shows such an example, which also confirms that the GC asymptotic solution for large k agrees well with the numerically obtained GC sampling structure from the time series. By the model construction (4), there is a causal flow from Y t to X t . However, through our asymptotic analysis, we confirm that the corresponding GC vanishes on certain sampling intervals, k.
To emphasize another kind of sampling artifact [30] in the GC sampling structure, we consider model II, for which β = + ∑ As shown in figure 1(d), model II also possesses properties I-II. Importantly, it displays a prominent feature of spurious GC as property x y k (the red curve in figure 1(d)) has a non-zero value, which becomes appreciable for > k 1, despite the fact that there is no causal influence from { } X t to { } Y t by the model construction.
This phenomenon is also discernible in figure 1(a) for the nonlinear network, for which can significantly deviate from zero over some interval τ, as clearly demonstrated by the 95% confidence region in figure 1(a).
As described above, properties I-III in the GC sampling structure are common features for both linear and nonlinear models. These properties have strong ramifications in the application of the GC analysis. Importantly, properties II and III demand additional sampling criteria for establishing the reliability of any conclusions about the causal influence for both linear and nonlinear systems.

Small τ limit for GC
Because GC is not invariant with respect to τ, we cannot choose the sampling-interval length arbitrarily in the GC analysis. One possible solution is to use a discrete time series sampled with very fine intervals from a continuous-time process. Intuitively, this approach would incorporate ever more information for causality determination with ever finer intervals. However, as shown in figure 2(a), for the case of nonlinear network dynamics, we have the following phenomenon as property (IV) the GC value, obtained using a discrete time series with the sampling interval τ, approaches zero almost linearly as τ → 0.
Therefore, the GC constructed in this way would give rise to a paradox that we eventually lose the ability for causal inference through GC despite the fact that more information is incorporated as τ → 0. We present the mathematical reasons below for property IV, from which one can see that it is also a general phenomenon.
Because GC can be computed through the spectral representation, we first study the limit of the spectral matrix ω τ ( )   have limits as τ → 0 and these limits are intrinsic properties of the continuous processes. Therefore, the GC is linearly proportional to the sampling interval length for small τ.  (6)). The inset is the coherence C (f) computed from PSD. Note that, by the asymptotic distribution theory of GC [27], the estimator of a directed GC has a bias p n / , where p is the regression order and n is the length of the time series. We have used the Bayesian information criterion [31] to determine the regression order p and have subtracted this type of bias in (a) and (b).

Procedure for reliable GC inference
From the above analysis, it is clear that one needs to be rather careful in interpreting causal inference using the GC analysis. Our analysis also provides a general approach to circumventing the issues discussed above: first, a range of τ should be used to sample continuous-time processes to obtain a set of discrete time series. Second, one computes the GC structure to examine its general behavior in order to (i) avoid using an accidental sampling interval τ for which GC may vanish despite the presence of causal influence, or (ii) avoid spurious non-zero GC values for finite τ as seen in figure 1(d) when there is no causal influence. Third, on account of the typical length scales τ osci and τ d of the oscillations and decay in the GC sampling structure, one should use fine sampling τ≪ τ osci , τ d to find the linear range of the GC, τ F ( ) (with the estimator bias [27] removed as in figure 2), as a function of τ, and then plot the ratio of τ τ F / ( ) to extract its limiting value as τ → 0. Returning to the network of two neurons, using the above approach, we can now successfully confirm the consistency between the causal connectivity and structural connectivity, i.e. that the topology = A 1 yx and = A 0 xy of the two-neuron networks can be as τ → 0, as shown in figure 2(b). In addition, we comment on some other aspects of the GC oscillation features observed in figure 1. As shown in figures 1(a) and (b) for both linear and nonlinear dynamics, the oscillation frequency in the GC sampling structure is twice the peak frequencies of the spectra. This relation also appears in the asymptotic solution (5) for model I, which contains the oscillation frequency β π = f / GC and whose ω ( ) S xy has a peak at . Note that for nonlinear networks, the GC analysis extracts an underlying regression structure. The spectra ( figure 1(a)) in the network case all have the same peak frequency, resembling those in model II ( figure 1(d)), therefore, they share similar GC sampling structures by the GC spectral theory. Complicated oscillation structures in the GC sampling structures may arise and, in general, are related to the spectral peak frequencies of the original time series.

An example: small-world networks
Finally, we illustrate the importance of the above GC extraction strategy in an example of the reconstruction of a directed, small-world network of neurons. The small world network is initially wired using a similar method to the Watts-Strogatz method [22] with a modification to take into account directed edges [32,33]. The small-world attributes are reflected in a large clustering coefficient and a small shortest path length in comparison with random networks. For our network of 100 excitatory neurons, the average clustering coefficient and the average shortest path lengths are 0.434 and 6.53, respectively. (For an undirected random network of the same number of nodes and edges, the average clustering coefficient and the average shortest path length would be 0.040, 2.44, respectively. For a directed random network, a node can then no longer always reach any other node through a directed path.) Using the voltage time series obtained from solving system (2), we can compute the conditional GC from neuron i to neuron j given the information about other neurons in the network (for details see [34]). The network topology is then successfully reconstructed (i.e. identical to that of the original network), as shown in figure 3(a), using the above GC extraction strategy. If a fixed sampling interval, say, τ = 6.0 ms, was used, the conventional GC analysis would produce a network topology that completely fails to capture the original network topology. This reconstructed network is shown in figure 3(b), which recovers only five directed connections while failing to find nearly 400 original connections in the total 9900 possible edges. Clearly, this reconstructed network, with a vanishing clustering coefficient and no connected paths, has even failed to capture the statistical structures of the network. Incidentally, we can also successfully reconstruct networks of coupled excitatory and inhibitory neurons with an arbitrary adjacent coupling matrix [19]. This example suggests that one should be very careful about the interpretation and implication of causal connectivity obtained using conventional methods for time series with a fixed sampling frequency. A systematic verification is needed to establish their validity.

Discussion and conclusion
In summary, we have shown that for both linear and nonlinear processes the computed GC is dependent on the sampling interval length τ. For instance, the GC may vanish on a certain set of τ despite the presence of true causal influence, or it can become non-zero over a range of τ despite the absence of any causal influence. Furthermore, the naive idea of using a sufficiently fine sampling interval, thus incorporating as much information as possible, will always give rise to vanishing GC regardless of whether there is a causal influence between the time series or not. By constructing simple idealized models, we have ascertained analytically the above phenomena as observed in both linear and nonlinear processes, and proposed approaches to  Figure 3. Loss of small-world attributes. Left panel: a small-world I&F network of 100 excitatory neurons whose directed connections are constructed using our GC extraction strategy with the smallest sampling interval τ = 1 ms (see text). These have been verified to be identical to the original directed network topology; here, a neuron is represented by a small circle and a directed edge is indicated by an arrow pointing from the presynaptic neuron to the postsynaptic neuron. Right panel: a network topology constructed using the τ = 6 ms sampling interval in the conventional GC analysis. This has recovered only five directed connections, while failing to capture nearly ∼400 original connections. This drastic loss of small-world attributes is a consequence of the incorrect sampling interval used in the GC analysis. The parameters for the I&F system (2) are μ = 1 kHz, λ = 0.012 and s = 0.005. overcome these sampling artifacts in order to obtain a reliable GC inference. We discussed the sampling issues of GC using the parametric methods of GC evaluations in the time domain as above. However, we point out that, naturally, similar issues about the sampling effects in the GC analysis also exist for non-parametric methods in the frequency domain. One has to deal carefully with the sampling issue regardless of whether a parametric method or non-parametric method of computing GC is used. We further comment that there is another related, interesting scientific question, namely, even for a fixed sampling rate, whether there is causal influence from one time scale to another in a coarse-grained sense. We note that this question was addressed in [35] by using a wavelet decomposition in time. Finally, we expect that our approach of obtaining reliable GC values can be used in analyzing GC relations for various kinds of observational data, such as EEG and fMRI signals, and may shed light on the impact of sampling effects in other empirical data-inference methods.