{'fraction_non_alphanumeric': 0.030123791977762694, 'fraction_numerical': 0.010472219528750121, 'mean_word_length': 4.555755072724008, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Granger Causality: Basic Theory and Application to Neuroscience

Multi-electrode neurophysiological recordings produce massive quantities of data. Multivariate time series analysis provides the basic framework for analyzing the patterns of neural interactions in these data. It has long been recognized that neural interactions are directional. Being able to assess the directionality of neuronal interactions is thus a highly desired capability for understanding the cooperative nature of neural computation. Research over the last few years has shown that Granger causality is a key technique to furnish this capability. The main goal of this article is to provide an expository introduction to the concept of Granger causality. Mathematical frameworks for both bivariate Granger causality and conditional Granger causality are developed in detail with particular emphasis on their spectral representations. The technique is demonstrated in numerical examples where the exact answers of causal influences are known. It is then applied to analyze multichannel local field potentials recorded from monkeys performing a visuomotor task. Our results are shown to be physiologically interpretable and yield new insights into the dynamical organization of large-scale oscillatory cortical networks.


Introduction
In neuroscience, as in many other fields of science and engineering, signals of interest are often collected in the form of multiple simultaneous time series. To evaluate the statistical interdependence among these signals, one calculates cross correlation functions in the time domain and ordinary coherence functions in the spectral domain. However, in many situations of interest, symmetric 1 measures like ordinary coherence are not completely satisfactory, and further dissection of the interaction patterns among the recorded signals is required to parcel out effective functional connectivity in complex networks. Recent work has begun to consider the causal influence one neural time series can exert on another. The basic idea can be traced back to Wiener [1] who conceived the notion that, if the prediction of one time series could be improved by incorporating the knowledge of a second one, then the second series is said to have a causal influence on the first. Wiener's idea lacks the machinery for practical implementation. Granger later formalized the prediction idea in the context of linear regression models [2]. Specifically, if the variance of the autoregressive prediction error of the first time series at the present time is reduced by inclusion of past measurements from the second time series, then the second time series is said to have a causal influence on the first one. The roles of the two time series can be reversed to address the question of causal influence in the opposite direction. From this definition it is clear that the flow of time plays a vital role in allowing inferences to be made about directional causal influences from time series data. The interaction discovered in this way may be reciprocal or it may be unidirectional.
Two additional developments of Granger's causality idea are important. First, for three or more simultaneous time series, the causal relation between any two of the series may be direct, may be mediated by a third one, or may be a combination of both. This situation can be addressed by the technique of conditional Granger causality. Second, natural time series, including ones from economics and neurobiology, contain oscillatory aspects in specific frequency bands. It is thus desirable to have a spectral representation of causal influence. Major progress in this direction has been made by Geweke [3,4] who found a novel time series decomposition technique that expresses the time domain Granger causality in terms of its frequency content. In this article we review the essential mathematical elements of Granger causality with special emphasis on its spectral decomposition. We then discuss practical issues concerning how to estimate such measures from time series data. Simulations are used to illustrate the theoretical concepts. Finally, we apply the technique to analyze the dynamics of a large-scale sensorimotor network in the cerebral cortex dur-ing cognitive performance. Our result demonstrates that, for a well designed experiment, a carefully executed causality analysis can reveal insights that are not possible with other techniques.

Bivariate Time Series and Pairwise Granger Causality
Our exposition in this and the next section follows closely that of Geweke [3,4]. To avoid excessive mathematical complexity we develop the analysis framework for two time series. The framework can be generalized to two sets of time series [3].

Time Domain Formulation
Consider two stochastic processes X t and Y t . Assume that they are jointly stationary. Individually, under fairly general conditions, each process admits an autoregressive representation (1) Jointly, they are represented as where the noise terms are uncorrelated over time and their contemporaneous covariance matrix is The entries are defined as Σ 2 = var(ǫ 2t ), Γ 2 = var(η 2t ), Υ 2 = cov(ǫ 2t , η 2t ). If X t and Y t are independent, then {b 2j } and {c 2j } are uniformly zero, Υ 2 = 0, Σ 1 = Σ 2 and Γ 1 = Γ 2 . This observation motivates the definition of total interdependence between X t and Y t as where | · | denotes the determinant of the enclosed matrix. According to this definition, F X,Y = 0 when the two time series are independent, and F X,Y > 0 when they are not.
Consider Eqs. (1) and (3). The value of Σ 1 measures the accuracy of the autoregressive prediction of X t based on its previous values, whereas the value of Σ 2 represents the accuracy of predicting the present value of X t based on the previous values of both X t and Y t . According to Wiener [1] and Granger [2], if Σ 2 is less than Σ 1 in some suitable statistical sense, then Y t is said to have a causal influence on X t . We quantify this causal influence by It is clear that F Y →X = 0 when there is no causal influence from Y to X and F Y →X > 0 when there is. Similarly, one can define causal influence from X to Y as It is possible that the interdependence between X t and Y t cannot be fully explained by their interactions. The remaining interdependence is captured by Υ 2 , the covariance between ǫ 2t and η 2t . This interdependence is referred to as instantaneous causality and is characterized by When Υ 2 is zero, F X·Y is also zero. When Υ 2 is not zero, F X·Y > 0.
The above definitions imply that Thus we decompose the total interdependence between two time series X t and Y t into three components: two directional causal influences due to their interaction patterns, and the instantaneous causality due to factors possibly exogenous to the (X, Y ) system (e.g. a common driving input).

Frequency Domain Formulation
To begin we define the lag operator L to be LX t = X t−1 . Rewrite Eqs. (3) and (4) in terms of the lag operator where a 2 (0) = 1, b 2 (0) = 0, c 2 (0) = 0, d 2 (0) = 1. Fourier transforming both sides of Eq. (11) leads to where the components of the coefficient matrix A(ω) are Recasting Eq. (12) into the transfer function format we obtain where the transfer function is H(ω) = A −1 (ω) whose components are After proper ensemble averaging we have the spectral matrix where * denotes complex conjugate and matrix transpose.
The spectral matrix contains cross spectra and auto spectra. If X t and Y t are independent, then the cross spectra are zero and |S(ω)| equals the product of two auto spectra. This observation motivates the spectral domain representation of total interdependence between X t and Y t as where |S(ω)| = S xx (ω)S yy (ω) − S xy (ω)S yx (ω) and S yx (ω) = S * xy (ω). It is easy to see that this decomposition of interdependence is related to coherence by the following relation: where coherence is defined as The coherence defined in this way is sometimes referred to as the squared coherence.
To obtain the frequency decomposition of the time domain causality defined in the previous section, we look at the auto spectrum of X t : It is instructive to consider the case where Υ 2 = 0. In this case there is no instantaneous causality and the interdependence between X t and Y t is entirely due to their interactions through the regression terms on the right hand sides of Eqs. (3) and (4). The spectrum has two terms. The first term, viewed as the intrinsic part, involves only the variance of ǫ 2t , which is the noise term that drives the X t time series. The second term, viewed as the causal part, involves only the variance of η 2t , which is the noise term that drives Y t . This power decomposition into an "intrinsic" term and a "causal" term will become important for defining a measure for spectral domain causality.
When Υ 2 is not zero it becomes harder to attribute the power of the X t series to different sources. Here we consider a transformation introduced by Geweke [3] that removes the cross term and makes the identification of an intrinsic power term and a causal power term possible. The procedure is called normalization and it consists of left-multiplying on both sides of Eq. (12). The result is where Since detÃ = det A we havẽ From the construction it is easy to see that E x andẼ y are uncorrelated, that is, cov(E x ,Ẽ y ) = 0. The variance of the noise term for the normalized Y t equation From Eq. (20), following the same steps that lead to Eq.
(18), the spectrum of X t is found to be: Here the first term is interpreted as the intrinsic power and the second term as the causal power of X t due to Y t . This is an important relation because it explicitly identifies that portion of the total power of X t at frequency ω that is contributed by Y t . Based on this interpretation we define the causal influence from Y t to X t at frequency ω as Note that this definition of causal influence is expressed in terms of the intrinsic power rather than the causal power. It is expressed in this way so that the causal influence is zero when the causal power is zero (i.e., the intrinsic power equals the total power), and increases as the causal power increases (i.e., the intrinsic power decreases).
By taking the transformation matrix as  and performing the same analysis, we get the causal influence from X t to Y t : By defining the spectral decomposition of instantaneous causality as [5] f X·Y (ω) = ln we achieve a spectral domain expression for the total interdependence that is analogous to Eq. (10) in the time domain, namely: We caution that the spectral instantaneous causality may become negative for some frequencies in certain situations and may not have a readily interpretable physical meaning.
It is important to note that, under general conditions, these spectral measures relate to the time domain measures as: The existence of these equalities gives credence to the spectral decomposition procedures described above.

Trivariate Time Series and Conditional Granger Causality
For three or more time series one can perform a pairwise analysis and thus reduce the problem to a bivariate problem. This approach has some inherent limitations. For example, for the two coupling schemes in Figure 1, a pairwise analysis will give the same patterns of connectivity like that in Figure 1(b). Another example involves three processes where one process drives the other two with differential time delays. A pairwise analysis would indicate a causal influence from the process that receives an early input to the process that receives a late input. To disambiguate these situations requires additional measures. Here we define conditional Granger causality which has the ability to resolve whether the interaction between two time series is direct or is mediated by another recorded time series and whether the causal influence is simply due to differential time delays in their respective driving inputs. Our development is for three time series. The framework can be generalized to three sets of time series [4].

Time Domain Formulation
Consider three stochastic processes X t , Y t and Z t . Suppose that a pairwise analysis reveals a causal influence from Y t to X t . To examine whether this influence has a direct component (Figure 1(b)) or is mediated entirely by Z t (Figure 1(a)) we carry out the following procedure. First, let the joint autoregressive representation of X t and Z t be where the covariance matrix of the noise terms is Next we consider the joint autoregressive representation of all three processes X t , Y t and Z t where the covariance matrix of the noise terms is From these two sets of equations we define the Granger causality from Y t to X t conditional on Z t to be The intuitive meaning of this definition is quite clear. When the causal influence from Y t to X t is entirely mediated by Z t (Fig. 1(a)), {b 4j } is uniformly zero, and Σ xx = Σ 3 . Thus, we have F Y →X|Z = 0, meaning that no further improvement in the prediction of X t can be expected by including past measurements of Y t . On the other hand, when there is still a direct component from Y t to X t ( Fig. 1(b)), the inclusion of past measurements of Y t in addition to that of X t and Z t results in better predictions of X t , leading to Σ xx < Σ 3 , and F Y →X|Z > 0.

Frequency Domain Formulation
To derive the spectral decomposition of the time domain conditional Granger causality we carry out a normalization procedure like that for the bivariate case. For Eqs. (29) and (30) the normalized equations are where (0) is generally not zero. We note that var(x * t ) = Σ 3 and this becomes useful in what follows.
For Eqs. (32), (33) and (34) the normalization process involves left-multiplying both sides by the matrix We denote the normalized equations as where the noise terms are independent, and their respective variances areΣ xx , Σ yy andΣ zz .
To proceed further we need the following important relation [4] and its frequency domain counterpart: To obtain f Y Z * →X * (ω), we need to decompose the spectrum of X * . The Fourier transform of Eqs. (36) and (37) gives: Assuming that X(ω) and Z(ω) from Eq. (40) can be equated with that from Eq. (41), we combine Eq. (40) and Eq. (41) to yield, where Q(ω) = G −1 (ω)H(ω). After suitable ensemble averaging, the spectral matrix can be obtained from which the power spectrum of X * is found to be The first term can be thought of as the intrinsic power and the remaining two terms as the combined causal influences from Y and Z * . This interpretation leads immediately to the definition We note that S x * x * (ω) is actually the variance of ǫ 3t as pointed out earlier. On the basis of the relation in Eq. (39), the final expression for Granger causality It can be shown that f Y →X|Z (ω) relates to the time domain measure F Y →X|Z via under general conditions. The above derivation is made possible by the key assumption that X(ω) and Z(ω) in Eq. (40) and in Eq. (41) are identical. This certainly holds true on purely theoretical grounds, and it may very well be true for simple mathematical systems. For actual physical data, however, this condition may be very hard to satisfy due to practical estimation errors. In a recent paper we developed a partition matrix technique to overcome this problem [6]. The subsequent calculations of conditional Granger causality are based on this partition matrix procedure.

Estimation of Autoregressive Models
The preceding theoretical development assumes that the time series can be well represented by autoregressive processes. Such theoretical autoregressive processes have infinite model orders. Here we discuss how to estimate autoregressive models from empirical time series data, with emphasis on the incorporation of multiple time series segments into the estimation procedure [7]. This consideration is motivated by the goal of applying autoregressive modeling in neuroscience. It is typical in behavioral and cognitive neuroscience experiments for the same event to be repeated on many successive trials. Under appropriate conditions, time series data recorded from these repeated trials may be viewed as realizations of a common underlying stochastic process.
Let X t = [X 1t , X 2t , · · · , X pt ] T be a p dimensional random process. Here T denotes matrix transposition. In multivariate neural data, p represents the total number of recording channels. Assume that the process X t is stationary and can be described by the following mth order autoregressive equation where A(i) are p × p coefficient matrices and E t = [E 1t , E 2t , · · · , E pt ] T is a zero mean uncorrelated noise vector with covariance matrix Σ.
To estimate A(i) and Σ, we multiply Eq. (46) from the right by X T t−k , where k = 1, 2, · · · , m. Taking expectations, we obtain the Yule-Walker equations where R(n) =< X t X T t+n > is X t 's covariance matrix of lag n. In deriving these equations, we have used the fact that < E t X T t−k >= 0 as a result of E t being an uncorrelated process.
For a single realization of the X process, {x i } N i=1 , we compute the covariance matrix in Eq. (47) according tõ If multiple realizations of the same process are available, then we compute the above quantity for each realization, and average across all the realizations to obtain the final estimate of the covariance matrix. (Note that for a single short trial of data one uses the divisor N for evaluating covariance to reduce inconsistency. Due to the availability of multiple trials in neural applications, we have used the divisor (N − n) in the above definition Eq. (48) to achieve an unbiased estimate.) It is quite clear that, for a single realization, if N is small, one will not get good estimates of R(n) and hence will not be able to obtain a good model. This problem can be overcome if a large number of realizations of the same process is available. In this case the length of each realization can be as short as the model order m plus 1.
Equations (46) contain a total of mp 2 unknown model coefficients. In (47) there is exactly the same number of simultaneous linear equations. One can simply solve these equations to obtain the model coefficients. An alternative approach is to use the Levinson, Wiggins, Robinson (LWR) algorithm, which is a more robust solution procedure based on the ideas of maximum entropy. This algorithm was implemented in the analysis of neural data described below. The noise covariance matrix Σ may be obtained as part of the LWR algorithm. Otherwise one may obtain Σ through Here we note that R T (k) = R(−k). where is the transfer function with A(0) being the identity matrix. From Eq. (52), after proper ensemble averaging, we obtain the spectral matrix Once we obtain the transfer function, the noise covariance, and the spectral matrix, we can then carry out causality analysis according to the procedures outlined in the previous sections.

Numerical Examples
In this section we consider three examples that illustrate various aspects of the general approach outlined earlier.

Example 1
Consider the following AR(2) model: where ǫ t , η t are gaussian white noise processes with zero means and variances σ 2 1 = 1, σ 2 2 = 0.7, respectively. The covariance between ǫ t and η t is 0.4. From the construction of the model, we can see that X t has a causal influence on Y t and that there is also instantaneous causality between X t and Y t .
We simulated Eq. (55) to generate a data set of 500 realizations of 100 time points each. Assuming no knowledge of Eq. (55) we fitted a MVAR model on the generated data set and calculated power, coherence and Granger causality spectra. The result is shown in Figure 2. The interdependence spectrum is computed according to Eq. (17) and the total causality is defined as the sum of directional causalities and the instantaneous causality. The result clearly recovers the pattern of connectivity in Eq. (55). It also illustrates that the interdependence spectrum, as computed according to Eq. (17), is almost identical to the total causality spectrum as defined on the right hand side of Eq. (27).

Example 2
Here we consider two models. The first consists of three time series simulating the case shown in Figure 1(a), in which the causal influence from Y t to X t is indirect and completely mediated by Z t : The second model creates a situation corresponding to Figure 1(b), containing both direct and indirect causal influences from Y t to X t . This is achieved by using the same system as in Eq. (56), but with an additional term in the first  Power (black for X, gray for Y ) spectra, interdependence spectrum (related to the coherence spectrum), and Granger causality spectra are displayed. Note that the total causality spectrum, representing the sum of directional causalities and the instantaneous causality, is nearly identical to the interdependence spectrum. equation: For both models. ǫ(t), ξ(t), η(t) are three independent gaussian white noise processes with zero means and variances of σ 2 1 = 0.3, σ 2 2 = 1, σ 2 3 = 0.2, respectively.
Each model was simulated to generate a data set of 500 realizations of 100 time points each. First, pairwise Granger causality analysis was performed on the simulated data set of each model. The results are shown in Figure 3(a), with the dashed curves showing the results for the first model and the solid curves for the second model. From these plots it is clear that pairwise analysis cannot differentiate the two coupling schemes. This problem occurs because the indirect causal influence from Y t to X t that depends completely on Z t in the first model cannot be clearly distinguished from the direct influence from Y t to X t in the second model. Next, conditional Granger causality analysis was performed on both simulated data sets. The Granger causality spectra from Y t to X t conditional on Z t are shown in Figure 3(b), with the second model's result shown as the solid curve and the first model's result as the dashed curve. Clearly, the causal influence from Y t to X t that was prominent in the pairwise analysis of the first model in Figure 3(a), is no longer present in Figure 3(b). Thus, by correctly determining that there is no direct causal influence from Y t to X t in the first model, the conditional Granger causality analysis provides an unambiguous dissociation of the coupling schemes represented by the two models.
We simulated the network model to generate a data set of 500 realizations each with 10 time points. Assuming no knowledge of the model, we fitted a 5th order MVAR model on the generated data set and performed power spectra, coherence and Granger causality analysis on the fitted model. The results of power spectra are given in the diagonal panels of Figure 4(b). It is clearly seen that all five oscillators have a spectral peak at around 25Hz and the fifth has some additional high frequency activity as well. The results of pairwise Granger causality spectra are shown in the off-diagonal panels of Figure 4(b) (solid curves). Compared to the network diagram in Figure 4(a) we can see that pairwise analysis yields connections that can be the result of direct causal influences (e.g. 1 → 2), indirect causal influences (e.g. 1 → 5) and differentially delayed driving inputs (e.g. 2 → 3). We further performed a conditional Granger causality analysis in which the direct causal influence between any two nodes are examined while the influences from the other three nodes are conditioned out. The results are shown as dashed curves in Figure  4(b). For many pairs the dashed curves and solid curves coincide (e.g. 1 → 2), indicating that the underlying causal influence is direct. For other pairs the dashed curves become zero, indicating that the causal influences in these pairs are either indirect are the result of differentially delayed inputs. These results demonstrate that conditional Granger causality furnishes a more precise network connectivity diagram that matches the known structural connectivity. One noteworthy feature about Figure 4(b) is that the spectral features (e.g. peak frequency) are consistent across both power and Granger causality spectra. This is important since it allows us to link local dynamics with that of the network.
Our results are the first to demonstrate in awake monkeys that synchronized beta oscillations not only bind multiple sensorimotor areas into a large-scale network during motor maintenance behavior, but also carry Granger causal influences from primary somatosensory and inferior posterior parietal cortices to motor cortex. Furthermore, the Granger causality graphs in Figure 5 provide a basis for fruitful speculation about the functional role of each cortical area in the sensorimotor network. First, steady pressure maintenance is akin to a closed loop control problem and as such, sensory feedback is expected to provide critical input needed for cortical assessment of the current state of the behavior. This notion is consistent with our observation that primary somatosensory area (S1) serves as the dominant source of causal influences to other areas in the network. Second, posterior parietal area 7b is known to be involved in nonvisually guided movement. As a higher-order association area it may maintain representations relating to the current goals of the motor system. This would imply that area 7b receives sensory updates from area S1 and outputs correctional signals to the motor cortex (M1). This conceptualization is consistent with the causality patterns in Figure 5. As mentioned earlier, previous work has identified beta range oscillations in the motor cortex as an important neural correlate of pressure maintenance behavior. The main contribution of our work is to demonstrate that the beta network exists on a much larger scale and that postcentral areas play a key role in organizing the dynamics of the cortical network. The latter conclusion is made possible by the directional information provided by Granger causality analysis. Since the above analysis was pairwise, it had the disadvantage of not distinguishing between direct and indirect causal influences. In particular, in monkey GE, the possibility existed that the causal influence from area S1 to inferior posterior parietal area 7a was actually mediated by inferior posterior parietal area 7b (Figure 5(a)). We used conditional Granger causality to test the hypothesis that the S1 → 7a influence was mediated by area 7b. In Figure 6(a) is presented the pairwise Granger causality spectrum from S1 to 7a (S1 → 7a, dark solid curve), showing significant causal influence in the beta frequency range. Superimposed in Figure 6(a) is the conditional Granger causality spectrum for the same pair, but with area 7b taken into account (S1 → 7a|7b, light solid curve). The corresponding 99% significance thresholds are also presented (light and dark dashed lines coincide). These significance thresholds were determined using a permutation procedure that involved creating 500 permutations of the local field potential data set by random rearrangement of the trial order independently for each channel (site). Since the test was performed separately for each frequency, a correction was necessary for the multiple comparisons over the whole range of frequencies. The Bonferroni correction could not be employed because these multiple comparisons were not independent. An alternative strategy was employed following Blair and Karniski [17]. The Granger causality spectrum was computed for each permutation, and then the maximum causality value over the frequency range was identified. After 500 permutation steps, a distribution of maximum causality values was created. Choosing a p-value at p = 0.01 for this distribution gave the thresholds shown in Figure 6 We see from Figure 6(a) that the conditional Granger causality is greatly reduced in the beta frequency range and no longer significant, meaning that the causal influence from S1 to 7a is most likely an indirect effect mediated by 7b. This conclusion is consistent with the known neuroanatomy of the sensorimotor cortex [18] in which area 7a receives direct projections from area 7b which in turn receives direct projections from the primary somatosensory cortex. No pathway is known to project directly from the primary somatosensory cortex to area 7a.
From Figure 5(a) we see that the possibility also existed that the causal influence from S1 to the primary motor cortex (M1) in monkey GE was mediated by area 7b. To test this possibility, the Granger causality spectrum from S1 to M1 (S1 → M1, dark solid curve in Figure 6(b)) was compared with the conditional Granger causality spectrum with 7b taken into account (S1 → M1|7b, light solid curve in Figure 6(b)). In contrast to Figure 6(a), we see that the beta-frequency conditional Granger causality in Figure 6(b) is only partially reduced, and remains well above the 99% significance level. From Figure 4(b), we see that the same possibility existed in monkey LU of the S1 to M1 causal influence being mediated by 7b. However, just as in Figure 6 is only partially reduced, and remains well above the 99% significance level.
The results from both monkeys thus indicate that the observed Granger causal influence from the primary somatosensory cortex to the primary motor cortex was not simply an indirect effect mediated by area 7b. However, we further found that area 7b did play a role in mediating the S1 to M1 causal influence in both monkeys. This was determined by comparing the means of bootstrap resampled distributions of the peak beta Granger causality values from the spectra of S1 → M1 and S1 → M1|7b by the Student's t-test. The significant reduction of beta-frequency Granger causality when area 7b is taken into account (t = 17.2 for GE; t = 18.2 for LU, p <<< 0.001 for both), indicates that the influence from the primary somatosensory to primary motor area was partially mediated by area 7b. Such an influence is consistent with the known neuroanatomy [18] where the primary somatosensory area projects directly to both the motor cortex and area 7b, and area 7b projects directly to primary motor cortex.

Summary
In this article we have introduced the mathematical formalism for estimating Granger causality in both time and spectral domain from time series data. Demonstrations of the technique's utilities are carried out on both simulated data, where the patterns of interactions are known, and on local field potential recordings from monkeys performing a cognitive task. For the latter we have stressed the physiological interpretability of the findings and pointed out the new insights afforded by these findings. It is our belief that Granger causality offers a new way of looking at cooperative neural computation and it enhances our ability to identify key brain structures underlying the organization of a given brain function.