Physiological Measurement
Related content
PAPER
Extended Granger causality: a new tool to identify the structure of physiological networks To cite this article: L Schiatti et al 2015 Physiol. Meas. 36 827
View the article online for updates and enhancements.
- Linear and non-linear brain–heart and brain–brain interactions during sleep L Faes, D Marinazzo, F Jurysta et al. - Cerebrovascular and cardiovascular variability interactions investigated through conditional joint transfer entropy in subjects prone to postural syncope Vlasta Bari, Beatrice De Maria, Claudio Enrico Mazzucco et al. - Basic cardiovascular variability signals: mutual directed interactions explored in the information domain Michal Javorka, Jana Krohova, Barbora Czippelova et al.
Recent citations - Inferring connectivity in networked dynamical systems: Challenges using Granger causality Bethany Lusch et al - Inferring Weighted Directed Association Network from Multivariate Time Series with a Synthetic Method of Partial Symbolic Transfer Entropy Spectrum and Granger Causality Yanzhu Hu et al - Inferring Weighted Directed Association Networks from Multivariate Time Series with the Small-Shuffle Symbolic Transfer Entropy Spectrum Method Yanzhu Hu et al
This content was downloaded from IP address 165.123.34.86 on 04/12/2017 at 15:56
Institute of Physics and Engineering in Medicine Physiol. Meas. 36 (2015) 827–843
Physiological Measurement doi:10.1088/0967-3334/36/4/827
Extended Granger causality: a new tool to identify the structure of physiological networks L Schiatti1, G Nollo1,3, G Rossato2 and L Faes1,3 1
BIOtech, Department of Industrial Engineering, University of Trento, Italy Department of Neurology, Sacro Cuore Hospital, Negrar, Verona, Italy 3 IRCS PAT-FBK, Trento, Italy 2
E-mail: [email protected] Received 3 October 2014, revised 23 December 2014 Accepted for publication 13 January 2015 Published 23 March 2015 Abstract
Granger causality (GC) is a very popular tool for assessing the presence of directional interactions between two time series of a multivariate data set. In its original formulation, GC does not account for zero-lag correlations possibly existing between the observed time series. In the present study we compare the GC with a novel measure, termed extended GC (eGC), able to capture instantaneous causal relationships. We present a two-step procedure for the practical estimation of eGC based on first detecting the existence of zero-lag correlations, and then assigning them to one of the two possible causal directions using pairwise measures of non-Gaussianity. The proposed method was validated in a simulation study, showing that the estimation procedure based on the extended representation overcomes the limits of the classic computation of GC, correctly detecting the presence and direction of zero-lag interactions and providing a meaningful causal interpretation based on the eGC. Then, GC and eGC were computed on the physiological variability series of heart period (HP), mean arterial pressure (AP) and cerebral blood flow velocity (FV) in ten subjects with postural related syncope (PRS), during different epochs of an head-up tilt test protocol. We found that both measures reflect the baroreflex impairment and the loss of cerebral autoregulation during pre-syncope. Furthermore, eGC analysis suggests that fast, within-beat effects between AP and FV variability contribute substantially to the mutual regulation of these physiological variables, and may play an important role in the impairment of cerebrovascular regulation associated with PRS.
0967-3334/15/040827+17$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
827
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Keywords: multivariate time series, instantaneous effects, cardiovascular interactions, cerebral autoregulation (Some figures may appear in colour only in the online journal) 1. Introduction The problem of inferring causality, i.e. directional interactions from a set of measured variables, has gained in the last years more and more importance in the study of physiological systems. One approach that has become increasingly popular is that introduced originally by Wiener (1956), and later formalized by Granger (1969) in terms of linear vector autoregressive (VAR) modeling of multivariate stochastic processes. According to Wiener–Granger causality (GC), a process X ‘Granger-causes’ another process Y in the presence of another relevant interacting (possibly multivariate) process Z, if, in an appropriate statistical sense, the past of X assists in predicting the present of Y beyond the degree to which the past of Y and Z already predicts the present of Y. It is important to remark that the identification of a Granger-causal interaction is not identical to the detection of a physically instantiated causal interactions in a system. Although the two descriptions are intimately related (Seth and Edelman 2007), physically instantiated causal structure can only be unambiguously identified by perturbing a system and observing the consequences (Pearl 2000). Nonetheless, GC is a pragmatic and well defined measure, and has delivered many insights into the functional connectivity of systems in a variety of fields, particularly in neuroscience and physiology (Bressler and Seth 2011, Schulz et al 2013). In its most traditional formulation, the measure of GC results from a VAR model representation of the observed processes, whereby the interactions between the past and the present of all processes are described in terms of linear regressions. This formulation presupposes that the considered model captures the entire interaction structure of the observed processes. If this is not the case, the causal interpretation provided by the regression coefficients may not be justified, and the GC may measure misleading patterns of causality. The classical strictly causal VAR model commonly adopted for GC estimation is limited in this respect, because its coefficients describe the time-lagged effects between the processes but do not account for the instantaneous (i.e. not lagged) effects. As a consequence of this incomplete description, the presence of significant instantaneous effects may lead the standard VAR models to provide wrong interpretations about causality. In fact, recent studies (Hyvarinen et al 2010, Faes 2014) demonstrated the adverse impact of neglecting instantaneous effects on causality analysis based on VAR modeling, providing theoretical arguments showing that the omission of zero-lag correlations from the model can change drastically the values of time-lagged coefficients, and thus of the causality measures based on the model. As a possible solution, these studies proposed the utilization of an extended model accounting for both instantaneous and lagged effects. The extended model was devised combining a classical VAR model and a model of the instantaneous causal interactions among the observed variables commonly known as structural equation model (SEM). Model estimation was performed using classic least-squares regression for the VAR part, and a more sophisticated method called LiNGAM (linear non-Gaussian acyclic model) (Shimizu et al 2006) which exploits non-Gaussianity to address the known identifiability problem of SEM models. The exploitation of combined VAR and SEM models for the study of a generalized form of GC which aggregates both instantaneous and time lagged-effects started with the works 828
L Schiatti et al
Physiol. Meas. 36 (2015) 827
of Hyvarinen et al (2008, 2010), and was employed later on to propose improved measures of causality in the frequency domain (Faes et al 2013a). In the present work we focus on the impact of neglecting instantaneous effects on the computation of the traditional time-domain measure of GC from X to Y, that is defined relating the prediction error variances of two linear regressions, respectively including and excluding the past of X in the prediction of the present of Y. This approach is ubiquitously used in time series analysis without considering the instantaneous effects among the modeled time series (Brovelli et al 2004, Wang et al 2007, Barnett and Seth 2014). Here we show that, since GC reflects the strictly causal VAR representation of the interactions, computing it when instantaneous effects are not trivial may lead to wrong conclusions about the true Granger-causal effects. We propose a measure of extended GC (eGC) that overcomes the limitations of GC by properly accounting for zerolag effects in the linear regressions. We also present a procedure for the practical estimation of eGC which addresses the problem of assigning instantaneous effects to the regression models. This procedure employs a two-step approach based on first estimating the existence of zero-lag correlations between two processes, and then orienting them along one of the two possible causal directions using pairwise measures of non-Gaussianity (Hyvarinen and Smith 2013). After validating the proposed method in a simulation study, we applied it to a practical case of study: the characterization of cardiovascular (CV) and cerebrovascular (CB) regulation mechanisms in subjects prone to postural-related syncope (PRS). PRS is an important and common clinical problem, consisting in a transient loss of consciousness and postural tone, with spontaneous recovery. Its underlying pathophysiological mechanisms are still not fully understood, although it has been associated with an impairment of CV and CB regulation (Toyry et al 1997, Dan et al 2002, Bechir et al 2003, Faes et al 2013a, 2013b). In this study we observe CV and CB variability measuring the fluctuations of heart period (HP), arterial pressure (AP), and cerebral blood flow velocity (FV) during a head-up tilt (HUT) testing protocol that is adopted in clinics to evoke PRS. Then, we compute traditional and novel measures of GC on the three time series measured during specific epochs of the test, in order to determine how the patterns of directional interaction between these variables are affected by the development of PRS. In this physiological application, the extension of the analysis to the inclusion of instantaneous effects is motivated by the fact that zero-lag correlations are expected to be relevant as they are related to fast (i.e. within-beat) mechanisms of CV and CB regulation. In CV regulation, within-beat interactions between systolic AP and HP are compatible with the short latencies, typically lower than 1 s, of the human cardiac baroreflex (Borst and Karemaker 1983); the presence of these interactions has been associated with the fast response of the parasympathetic component of the autonomic nervous system, and has been studied during postural stress and PRS (Westerhof et al 2006, Faes et al 2013b). In CB regulation, studies of dynamic cerebral autoregulation have shown that the response of cerebral blood FV to fast transient changes in AP is completed within a few seconds, and occurs in large part before the next cardiac pulsation (Zhang et al 1998, Panerai et al 1999). The paper is structured as follows: in section 2 the proposed framework for the eGC estimation based on VAR modeling is presented; in section 3 a simulation study to test the effectiveness of eGC measure when instantaneous effects are present is shown; in section 4 the application of the proposed method to a real data set is described; section 5 concludes the paper. We distribute a Matlab toolbox for the computation of GC and eGC; the toolbox is available online at www.lucafaes.net/eGC.html.
829
L Schiatti et al
Physiol. Meas. 36 (2015) 827
2. Methods 2.1. Extended GC
GC is assessed in its traditional formulation based on VAR models (Bressler and Seth 2011). Specifically, given a multivariate stochastic process Y composed of M scalar processes Y1,...,YM, their dynamical interactions are described by the VAR model of order p p
Y (n ) =
∑ A (k ) Y (n − k )
+ U (n ) ,
(1)
k=1
where Y(n) = [Y1(n)⋅⋅⋅YM(n)]T is the vector variable obtained sampling the process Y at the time instant n, A(k) are M × M coefficient matrices in which the element Aji(k) describes the dependence of Yj(n) on Yi(n − k) (i,j = 1,...M; k = 1,...,p), and U = [U1⋅⋅⋅UM]T represents the model residuals. To compute GC from the ith process Yi to the jth process Yj, two linear regressions are considered starting from the VAR representation in (1). First, Yj is described in terms of the past of all processes in Y by means of a so-called unrestricted regression that is obtained considering the part of the joint representation in (1) relevant to the jth process p
Y j (n ) =
∑ Aj (k ) Y (n − k )
+ Uj (n ) ,
(2)
k=1
where Aj(k) is the jth row of A(k). Then, Yj is described in terms of the past of all processes in Y except Yi by means of the so-called restricted regression p
Y j (n ) =
∼
∑ Aj (k ) Y i (n − k )
∼ + Uj (n ) ,
(3)
k=1
where Yi = [Y1⋅⋅⋅Yi–1Yi+1⋅⋅⋅YM]T denotes the original vector process devoid of the ith scalar ∼ ∼ process. Note that the coefficients Aj and the residuals Uj yielded by the restricted regression (3) are in general different than the Aj and the Uj yielded by the unrestricted regression (2). The GC measure is computed by comparing the variance of the residuals resulting from the two regressions: ∼ var Uj (n ) GC . i → j = ln (4) var Uj (n ) The GC measure (4) is always positive and proportional to the extent to which the knowledge of the past of Yi helps improving the prediction of Yj(n) above and beyond the improvement achieved using the past of Yi; the measure becomes null if the past of Yi is ineffective to the prediction improvement (in this case, the coefficients Aji(k) are zero for all lags k and the restricted and unrestricted regressions (2) and (3) coincide). An important property of the GC is that it is sensitive only to the causal effects from Yi to Yj which are not mediated by the interaction of Yi and Yj with any other process Yk (k ≠ i, k ≠ j). Therefore, the GC is able to elicit direct dependencies between two processes and to discard the confounding effects of third processes belonging to the set of processes used for VAR modeling. The model of equation (1) is a strictly causal model, in the sense that it describes only the time-lagged interactions between the processes. Therefore, GC is computed without considering the instantaneous (i.e. non-delayed) effects among the observed time series. However, the possible presence of zero-lag effects is known to have an impact also on the time-lagged effects (Hyvarinen et al 2010, Faes 2014), and thus may affect the reliability of the observed 830
L Schiatti et al
Physiol. Meas. 36 (2015) 827
GC patterns. To overcome this problem, we propose the utilization of an extended version of GC, derived from an extended VAR model combining both instantaneous and lagged effects: p
Y (n ) =
∑ B (k ) Y (n − k )
+ W (n ) ,
(5)
k=0
where B(k) are M × M coefficient matrices in which the element Bji(k) describes the dependence of Yj(n) on Yi(n − k) (i,j = 1,...M; k = 0,...,p), including also coefficients weighing the present samples of Y, and W = [W1⋅⋅⋅WM]T represents the model residuals. To compute eGC from Yi to Yj, unrestricted and restricted linear regressions incorporating zero-lag causal effects are computed starting from the extended VAR representation in (5), analogously to what seen for GC: p
Y j (n ) =
∑ Bj (k ) Y (n − k )
+ Wj (n ) ,
(6)
+ W͠ j (n ) ,
(7)
k=0 p
Y j (n ) =
∼
∑ Bj (k ) Y i (n − k ) k=0
∼ where Bj(k) and Bj (k ) are row coefficients weighting respectively the present and past samples of all processes in Y for the unrestricted regression (6), and the present and past samples of all processes in Y except Yi for the restricted regression (7). The eGC measure is computed by comparing the variance of the residuals of the two regressions: eGC i → j = ln
var W͠ j (n ) var Wj (n )
.
(8)
2.2. Estimation
In practical analysis, GC and eGC are estimated between any pair of time series taken from an observed multivariate dataset by identifying the unrestricted and restricted regression models (respectively, equations (2) and (6) and equations (3) and (7)) and then combining the variance of the estimated model residuals as in equations (4) and (8). In this study, regressions were performed by the standard least squares approach, and the number of lags p was set identical for all time series by applying the Bayesian Information Criterion to the overall VAR representation of equation (1) (Faes et al 2012). Moreover, the statistical significance of each computed GC and eGC measure was assessed using the Fisher F-test. While for the computation of the traditional GC the estimation is straightforward as just described above, for the eGC analysis the procedure is complicated by the need of identifying the structure of the zero-lag effects among the time series. Such a structure cannot be imposed arbitrarily, because it has to be consistent with the joint representation of the observed multivariate time series (Shimizu et al 2006). The zero-lag effects are collected in the matrix B(0) of the extended VAR representation in (5): allowing an effect at lag zero between two processes Yi and Yj corresponds to set either Bij(0) ≠ 0 or Bji(0) ≠ 0; the coefficients cannot be both nonzero because the presence of a pairwise loop would make the extended model unidentifiable (Shimizu et al 2006, Hyvarinen et al 2010, Faes et al 2013a). Thus, the estimation of eGC requires to identify the presence and causal direction of the zero-lag effects between the observed processes to be incorporated into the regression models (6) and (7). In the following, we propose a two-step procedure to accomplish this task, based first on estimating the 831
L Schiatti et al
Physiol. Meas. 36 (2015) 827
existence of zero-lag correlations in an undirected sense, and then on finding their directions using pairwise measures of non-Gaussianity (Hyvarinen and Smith 2013). As a first step, we determine whether an effect at lag zero is present between any pair of processes Yi and Yj, looking at the correlation between the variables Yi(n) and Yj(n). The simplest way to assess correlations in a multivariate process is to compute its partial correlation matrix, which is defined as the inverse of the covariance matrix of the process (Marrelec et al 2006). Accordingly, in our preliminary study (Schiatti et al 2014) we computed the partial correlation of Y, [covY(n)]–1. However, the partial correlation reflects correctly the true correlations in the data only in the absence of any temporal structure, while it may lead to detect spurious instantaneous correlations when lagged correlations are also present. This occurs due to possible common driver effects from some past lag component, say Yk(n − l), which drives both Yi(n) and Yj(n) inducing a zero-lag correlation that is not conditioned out by the partial correlation. On the other hand, assuming that the VAR model (1) properly captures all time-lagged effects into its coefficients, the residuals U do not have any temporal structure, so that finding nonzero partial correlation between Ui(n) and Uj(n) means that a ‘true’ direct effect exists between Yi(n) and Yj(n). Therefore in this study we assessed zero-lag correlations in our multivariate process computing the partial correlation of the VAR model residuals, [covU(n)]−1. In the practical computation, significant nonzero values of the estimated partial correlations were determined using bootstrap resampling (Zoubir and Boashash 1998). Specifically, we drew 100 bootstrap samples of the same length from the original estimates of U(n) through sampling with replacement and computed the partial correlation for each sample. We then exploited the resulting distribution of values to get a 95% confidence interval for [covU(n)]−1, and labeled as statistically significant all elements of the original partial correlation matrix for which the confidence interval did not encompass the zero level. In the second step, the direction of the zero-lag interaction between pairs of series associated with residuals showing a significant partial correlation was assessed employing the pairwise likelihood ratios for non-Gaussian data recently proposed by Hyvärinen and Smith (2013). This approach is essentially based on estimating the likelihood of the models corresponding to the two directions of causal influence between a pair of variables, and on deciding in favor of the model with higher likelihood after the computation of a likelihood ratio. The two models, which consider the standardized (i.e. having zero mean and unit variance) variables x and y, are defined by the simple regressions y = ρx + d (model x → y) and x = ρy + e (model y → x), where d and e are the residuals, assumed to be independent of the regressors x and y, and ρ is the correlation coefficient between x and y, and is therefore the same in the two models. Considering N realizations of x and y, the ratio between the likelihoods of the two models, L(x → y) and L(y → x), is first computed as R =
L (x → y ) 1 ; log L (y → x ) N
(9)
then, the causal direction is chosen as x → y if R is positive, and as y → x if R is negative. In (Hyvarinen and Smith 2013), exploiting the information theoretic interpretation of the likelihood ratio and the connection between information measures and non-Gaussianity, it is shown that choosing the model with higher likelihood corresponds to choosing the direction which leads to maximum non-Gaussianity. This is in line with the knowledge that the direction of SEM models describing the instantaneous interactions among processes can be inferred only if the observed variables have a non-Gaussian distribution (Shimizu et al 2006, Hyvarinen et al 2010, Faes et al 2013a). For the practical computation of the likelihood ratio, we adopted the first-order approximation of the likelihood of the two models, describing the 832
L Schiatti et al
Physiol. Meas. 36 (2015) 827
log-probability density functions with the logistic density (Hyvarinen and Smith 2013). This leads to estimating the likelihood ratio as R ̂ = ρ ̂ x tanh(y )− y tanh(x ) ,
(10)
where ρ ̂ is the estimated correlation coefficient and denotes sample average. In our application, the likelihood ratio was computed as in (10) taking as variables x and y the standardized version of any pair of estimated VAR residuals Ui(n) and Uj(n) for which a statistically significant partial correlation was found according to bootstrap resampling. 3. Simulation study 3.1. Theoretical formulation
In order to describe the features of GC and eGC, and to illustrate their relationship, we consider the VAR process of order p = 2, composed by M = 3 scalar processes, generated by the equations:
⎧Y1(n )= δY2(n ) + b11(1)Y1(n − 1) + b11(2)Y1(n − 2) + W1(n ) ⎪ ⎨Y2(n )= b 22(1) Y2(n − 1) + b 22(2) Y2(n − 2) + b 21(2)Y1(n − 2) + b 23(2)Y3(n − 2) + W2(n ) , (11) ⎪Y (n )= δY (n ) + b (1)Y (n − 1) + W (n ) ⎩ 3 1 31 1 3
where the Wi are white and uncorrelated model residuals. The coefficients appearing on the right-hand side of the equations can be collected in the matrices B(k), k = 0,1,2. The diagonal elements B(k) are set in order to generate complex conjugate poles for the processes Y1 and Y2 (b11(1) = 2r1cos(2πf1); b11 (2) = −r12; b22(1) = 2r2cos(2πf2); b 22 (2) = −r22); we set modulus r1 = 0.9 and r2 = 0.8, and phase ±3/5π and ±π/5 determining autonomous oscillations at f1 = 0.3 Hz and f2 = 0.1 Hz, respectively for Y1 and Y2. Then, the off-diagonal elements of B(k) determine direct causal effects from one process to another. The causal effects at lag zero are set according to the parameter δ; here we consider the conditions δ = 0 and δ = 0.8, entailing absence and presence of instantaneous effects, respectively. The process generated by (11) is an extended VAR process, which can be described in a compact form as in equation (5). The corresponding strictly causal VAR representation can be obtained analytically from the extended representation using the relations (11): U (n ) = [I − B(0)]−1 W (n ),
(12)
A(k ) = [I − B(0) ]−1 B(k ) for each k ≥ 1.
(13)
These relations allow to compute the model coefficients and residuals of the strictly causal VAR model equivalent to the extended model defined by (11). Doing so, we depicted the patterns of causality resulting from the two representations in the two conditions δ = 0 and δ = 0.8, obtaining the connectivity diagrams of figure 1. When δ = 0 all the imposed effects from one process to another are lagged, and result from the nonzero model coefficients: Y1 → Y2 (b21(2) = 0.5), Y1 → Y3 (b31(1) = 0.5) and Y3 → Y2 (b23(2) = 0.5). In this condition with absence of instantaneous causality the extended and strictly causal VAR representations coincide (B(0) = 0 in (12) and (13) makes U(n) = W(n) and A(k) = B(k)), so that both models provide a correct interpretation of the causality patterns (figure 1(a)). When δ = 0.8, two instantaneous effects over the directions Y2 → Y1 and Y1 → Y3 are imposed in addition to the lagged causal effects described so far. In this case the extended VAR representation correctly reflects relations imposed by equation (11), as depicted in figure 1(b) (left). On the contrary, due to the 833
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Figure 1. Extended VAR representation (left) and strictly causal VAR representation (right) of the VAR processes generated by the equations in (11), in absence (a) and presence (b) of instantaneous effects. Arrows represent causal effects, with the numbers indicating the lag at which effects are detected.
presence of instantaneous causality the causal interpretation based on the strictly causal model is now different, as it reflects the strictly causal model coefficients computed using (13). In particular, we see that the strictly causal representation identifies more lagged relations than the imposed ones, interpreting as lagged the instantaneous effect Y2 → Y1 and detecting the spurious lagged causality effects Y3 → Y1 and Y2 → Y3, as well as an unexpected effect from Y3(n − 2) to Y3(n). This example shows that the traditional strictly causal representation may provide misleading patterns of causality when the instantaneous effects among the observed processes are not negligible. 3.2. Implementation on simulations
The feasibility of the estimation approaches illustrated in section 2.2 for the computation of GC and eGC measures was tested on simulations of the extended VAR process illustrated in section 3.1. Specifically, realizations of (11) lasting N = 300 points were generated starting from realizations of the extended residuals Wi, and the comparative analysis was performed under three different scenarios: (a) absence of instantaneous effects (δ = 0 in (11)); (b) presence of instantaneous effects (δ = 0.8) with non-Gaussian extended residuals; and (c) presence of instantaneous effects (δ = 0.8) with Gaussian extended residuals. The extended residuals, Wm(n), m = 1,…,M, n = 1,…,N, were generated as independent Gaussian white noises for scenarios (a) and (c), and as independent non-Gaussian white noises for scenario (b); in this last case, non-Gaussianity was achieved first generating Zm(n) as independent Gaussian white noises and then applying the nonlinear transformation Wm(n) = sign(Zm(n))|Zm(n)|q, with exponent q chosen in the range [0.5, 0.8] or [1.2, 2.0] to yield respectively a sub-Gaussian or super-Gaussian distribution for Wi (Faes et al 2013a). GC and eGC were computed from 100 realizations of the simulation generated as described above for the three scenarios. Moreover, the theoretical values of GC and eGC were also computed exploiting the procedure described in (Faes et al 2015, Barnett and Seth 2014). This procedure allows to express the variance of the residuals of both unrestricted and restricted models, and thus to express the GC and eGC, in terms of the VAR parameters. Here, theoretical 834
L Schiatti et al
Physiol. Meas. 36 (2015) 827
values of GC and eGC were used to assess the reliability of the measures estimated from the process realizations. Figure 2 reports the distributions of GC and eGC estimated over 100 realizations of the simulation under the three scenarios described above. Moreover, table 1 reports the number of realizations for which instantaneous effects were detected along any specified causal direction. In the absence of instantaneous effects, the distributions of the estimated GC and eGC overlap almost perfectly, and are centered in the theoretical values (figure 2(a)). Moreover, both estimates were always detected as statistically significant, according to the F-test, when computed along the coupled directions Y1 → Y2, Y1 → Y3 and Y3 → Y2, and were almost never detected as significant when computed along the uncoupled directions. This demonstrates that both approaches are able to retrieve the patterns of GC when the investigated processes do not interact at lag zero. The equivalence between GC and eGC in this case is confirmed by the observation that the bootstrap procedure performed to test for instantaneous causality detected the presence of significant partial correlations among the VAR residuals in a very limited number of realizations (table 1(a)). When significant instantaneous effects were simulated imposing δ = 0.8 in (11), the estimates of GC and eGC were still centered on their theoretical values (figure 2(b)), documenting the appropriateness of the estimation procedures. However, in this case the patterns of GC and eGC were substantially different: eGC was markedly higher than GC along the directions Y2 → Y1 and Y1 → Y3, and markedly lower and non-significant along the directions Y3 → Y1 and Y2 → Y3. Comparing figure 2(b) with figure 1(b), we see that the two measures reflect the extended and strictly causal representations of the simulated processes: the eGC is significant whenever an instantaneous or time lagged effect is set from one process to another, while the GC follows the misleading patterns of lagged causality provided by the strictly causal representation. This confirms the unsuitability of GC to measure Granger-causal effects in processes exhibiting zero-lag correlations. On the contrary, the estimation procedure based on the extended VAR representation correctly detects the presence and direction of zero-lag interactions among the processes (table 1(b)) and provides a meaningful causal interpretation based on the eGC (figure 2(b)). Figure 2(c) illustrates what happens when the requirement of non-Gaussianity of model residuals is not met. In this case, instantaneous effects can be detected if not present, or be assigned to the wrong causal direction (table 1(c)). As a consequence the distributions of the estimated eGC became wider and deviated from the theoretical values (figure 2(c)). This documents the importance of satisfying the requirement of non-Gaussianity to guarantee a proper identification of instantaneous effects, and consequently a reliable estimation of Grangercausal effects. 4. Real data application This section describes the application of the proposed framework based on VAR modeling to estimate patterns of GC and eGC in physiological time series measured in the context of CV and CB regulation in subjects with PRS. 4.1. Experimental protocol
Ten young subjects (3 males, 22.4 ± 9.5 years) with history of recurrent unexplained syncope were studied during a HUT test protocol. The acquired signals were the electrocardiogram (ECG), the finger photoplethysmographic AP, and the cerebral blood flow based on 835
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Figure 2. Box-plot of GC and eGC values computed over 100 realizations of equation (11) generated according to the three scenarios: (a) δ = 0, Gaussian Wi; (b) δ = 0.8, non-Gaussian Wi; (c) δ = 0.8, Gaussian Wi. The symbols ✠ indicate the theoretical value of GC and eGC. The number of GC or eGC values detected as statistically significant (F-test with 0.01 significance) is reported below each distribution.
transcranial Doppler ultrasonography, measured with subjects lying in the resting supine position and standing in the 60° upright position reached after passive HUT. Postural stress in the upright position was prolonged until the occurrence of PRS. All subjects exhibited spontaneous recovery with return to the supine position. The analyzed time series were the beat-to-beat variability of the HP, the mean AP and the mean cerebral blood FV, measured respectively as the temporal distance between two consecutive R peaks of the ECG, and as the integral of the AP and cerebral flow signals taken within each detected diastolic pulse interval normalized to interval duration (figure 3). Three stationary epochs of N = 300 beats were considered for the analysis, one in the supine position (SU), and two in the early (ET) and late (LT) phases of passive standing in the upright position (respectively at ~2 min after tilting and immediately before the drop in AP indicating presyncope). Details about experimental protocol can be found in (Faes et al 2013c). 836
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Table 1. Percentage of realizations of equation (11) for which instantaneous correlations
were detected along any given causal direction, looking at the statistical significance of the partial correlation of the VAR residuals, under the three scenarios: (a) δ = 0, Gaussian Wi; (b) δ = 0.8, non-Gaussian Wi; (c) δ = 0.8, Gaussian Wi. Zero-lag effects
(a) (b) (c)
Y2 → Y1
Y1 → Y2
Y3 → Y1
Y1 → Y3
Y3 → Y2
Y2 → Y3
6 100 52
1 0 48
1 0 46
2 100 54
3 0 4
3 8 5
Figure 3. HP, AP and FV variability series measure; (a) computation of HP from ECG signal; (b) computation of AP from finger photoplethysmographic AP signal; (c) computation of FV from transcranial Doppler ultrasonography signal.
4.2. Data analysis
GC and eGC were computed between each pair of the three considered time series in the SU, ET and LT windows. The model order was assessed using the Bayesian information criterion (Faes et al 2012), minimizing the corresponding figure of merit in the range 1–20. The identification of zero-lag correlations was performed using 100 bootstrap samples and a 95th percentile confidence interval. The F-test for assessing nonzero values of GC and eGC was run with 1% significance. The statistical significance of the differences in the values of GC or eGC across the three conditions (SU, ET, LT) was assessed by the one-way ANOVA with 5% significance, followed by post-hoc comparison based on the Student t-test for paired data with Bonferroni correction. Moreover, the statistical significance of the difference between the distributions of GC and eGC computed for a specific condition was assessed by the Student t-test for paired data. In addition, since the identification of instantaneous effects relies on the hypothesis of nonGaussianity for the model residuals, such an hypothesis was tested by means of a method based on mutual information (MI) (Hlinka et al 2011). This method exploits the knowledge that the MI I(x,y) between two standardized variables x and y has a lower bound such that I(x,y) ≥ IGauss(r) 837
L Schiatti et al
Physiol. Meas. 36 (2015) 827
= −0.5ln(1−r2), where r is the linear correlation between x and y; the equality holds exactly for a bivariate Gaussian distribution, for which the MI IGauss is a function of the correlation between the two variables. This allows to assess the deviation from Gaussianity of the joint distribution between x and y by comparing the MI I(x,y) with the MI between surrogates of x and y which have the same linear correlation of the original variables but are enforced to have bivariate Gaussian distribution. In our study, the role of x and y was taken by any pair of estimated VAR residuals Ui(n) and Uj(n), and surrogates were constructed as multivariate Fourier transform (FT) surrogates (Prichard and Theiler 1994). These surrogates were obtained computing the FT of the original series, keeping the magnitudes of the Fourier coefficients but adding the same random number to their phase, and finally returning into the time domain through inverse FT computation. We generated 100 surrogate pairs for each original pair of residuals, and rejected the null hypothesis of two linear Gaussian processes when the original MI was larger than the 95th percentile of the distribution of MI computed over the surrogate pairs. The MI was computed using the nearest neighbor estimator proposed by Kraskov et al (2004). 4.3. Results and discussion
The order of the regression models employed for the computation of GC and eGC, selected with the Bayesian Information Criterion, was comprised for all subjects between 2 and 5 and did not vary significantly across the three considered epochs of HUT (supine: 3.6 ± 1.1; early tilt: 3.1 ± 0,6; late tilt: 3.4 ± 0.5 mean ± SD across subjects). The analysis of GC and eGC, performed in terms of both magnitude of the measure and number of subjects with significant causality, led to the results depicted in figure 4. This analysis indicated the existence of a bidirectional causality between HP and AP, with prevalence of HP → AP causality, and of mostly unidirectional causality along the directions HP → FV and AP → FV. Moreover we found significantly higher causality AP → HP during ET compared with SU, and higher causality AP → FV during LT compared with SU. These results are in line with expected mechanisms of CV and CB regulation. Indeed, the bidirectional HP–AP interaction reflects the closed loop regulation of the CV system, resulting from the coexistence between the feedback baroreflex mechanism acting from AP to HP and the feedforward mechanical effects of HP on AP (Cohen and Taylor 2002). The strengthening of the causal coupling from AP to HP with the tilt transition documents the well known activation of the baroreflex (Porta et al 2011, Faes et al 2013b), developed as a reaction to the change of posture finalized to maintain blood pressure to a constant level; on the opposite we observed a decrease in the baroreflex coupling AP → HP moving from ET to LT, though not statistically significant, that might reflect an impairment of this mechanism in proximity of syncope (Faes et al 2006). As to the analysis of CB regulation, the significant causal effect of the two CV variables on FV are expected since it is known that the cerebral blood flow can vary rapidly in response to systemic AP variations (Zhang et al 2000). On the other hand, it is also known that the CB bed works for maintaining a relatively constant blood flow in the face of the perturbing action of AP, according to a homeostatic mechanism known as cerebral autoregulation (Paulson et al 1990). According to this mechanism, increases in the coupling between AP and FV have been often associated with a loss of cerebral autoregulation (Panerai et al 1998, Mitsis et al 2004), and in this view the higher GC from AP to FV which we observed upon presyncope may be indicative of impaired CB control. The results described above were confirmed using eGC, which indeed increased from SU to ET in the direction AP → HP (though without reaching statistical significance, p = 0.114), and increased significantly from SU to LT in the direction AP → FV. Furthermore, eGC analysis showed modifications of the CB causality patterns associated with presyncope that were 838
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Figure 4. Box-plot of GC and eGC values between HP, mean AP and mean FV series computed during SU, ET and LT epochs of HUT. The number of GC or eGC values detected as statistically significant (F-test with 0.01 significance) is reported below each distribution. Significant difference between pairs of distributions: *p
Related content
PAPER
Extended Granger causality: a new tool to identify the structure of physiological networks To cite this article: L Schiatti et al 2015 Physiol. Meas. 36 827
View the article online for updates and enhancements.
- Linear and non-linear brain–heart and brain–brain interactions during sleep L Faes, D Marinazzo, F Jurysta et al. - Cerebrovascular and cardiovascular variability interactions investigated through conditional joint transfer entropy in subjects prone to postural syncope Vlasta Bari, Beatrice De Maria, Claudio Enrico Mazzucco et al. - Basic cardiovascular variability signals: mutual directed interactions explored in the information domain Michal Javorka, Jana Krohova, Barbora Czippelova et al.
Recent citations - Inferring connectivity in networked dynamical systems: Challenges using Granger causality Bethany Lusch et al - Inferring Weighted Directed Association Network from Multivariate Time Series with a Synthetic Method of Partial Symbolic Transfer Entropy Spectrum and Granger Causality Yanzhu Hu et al - Inferring Weighted Directed Association Networks from Multivariate Time Series with the Small-Shuffle Symbolic Transfer Entropy Spectrum Method Yanzhu Hu et al
This content was downloaded from IP address 165.123.34.86 on 04/12/2017 at 15:56
Institute of Physics and Engineering in Medicine Physiol. Meas. 36 (2015) 827–843
Physiological Measurement doi:10.1088/0967-3334/36/4/827
Extended Granger causality: a new tool to identify the structure of physiological networks L Schiatti1, G Nollo1,3, G Rossato2 and L Faes1,3 1
BIOtech, Department of Industrial Engineering, University of Trento, Italy Department of Neurology, Sacro Cuore Hospital, Negrar, Verona, Italy 3 IRCS PAT-FBK, Trento, Italy 2
E-mail: [email protected] Received 3 October 2014, revised 23 December 2014 Accepted for publication 13 January 2015 Published 23 March 2015 Abstract
Granger causality (GC) is a very popular tool for assessing the presence of directional interactions between two time series of a multivariate data set. In its original formulation, GC does not account for zero-lag correlations possibly existing between the observed time series. In the present study we compare the GC with a novel measure, termed extended GC (eGC), able to capture instantaneous causal relationships. We present a two-step procedure for the practical estimation of eGC based on first detecting the existence of zero-lag correlations, and then assigning them to one of the two possible causal directions using pairwise measures of non-Gaussianity. The proposed method was validated in a simulation study, showing that the estimation procedure based on the extended representation overcomes the limits of the classic computation of GC, correctly detecting the presence and direction of zero-lag interactions and providing a meaningful causal interpretation based on the eGC. Then, GC and eGC were computed on the physiological variability series of heart period (HP), mean arterial pressure (AP) and cerebral blood flow velocity (FV) in ten subjects with postural related syncope (PRS), during different epochs of an head-up tilt test protocol. We found that both measures reflect the baroreflex impairment and the loss of cerebral autoregulation during pre-syncope. Furthermore, eGC analysis suggests that fast, within-beat effects between AP and FV variability contribute substantially to the mutual regulation of these physiological variables, and may play an important role in the impairment of cerebrovascular regulation associated with PRS.
0967-3334/15/040827+17$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
827
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Keywords: multivariate time series, instantaneous effects, cardiovascular interactions, cerebral autoregulation (Some figures may appear in colour only in the online journal) 1. Introduction The problem of inferring causality, i.e. directional interactions from a set of measured variables, has gained in the last years more and more importance in the study of physiological systems. One approach that has become increasingly popular is that introduced originally by Wiener (1956), and later formalized by Granger (1969) in terms of linear vector autoregressive (VAR) modeling of multivariate stochastic processes. According to Wiener–Granger causality (GC), a process X ‘Granger-causes’ another process Y in the presence of another relevant interacting (possibly multivariate) process Z, if, in an appropriate statistical sense, the past of X assists in predicting the present of Y beyond the degree to which the past of Y and Z already predicts the present of Y. It is important to remark that the identification of a Granger-causal interaction is not identical to the detection of a physically instantiated causal interactions in a system. Although the two descriptions are intimately related (Seth and Edelman 2007), physically instantiated causal structure can only be unambiguously identified by perturbing a system and observing the consequences (Pearl 2000). Nonetheless, GC is a pragmatic and well defined measure, and has delivered many insights into the functional connectivity of systems in a variety of fields, particularly in neuroscience and physiology (Bressler and Seth 2011, Schulz et al 2013). In its most traditional formulation, the measure of GC results from a VAR model representation of the observed processes, whereby the interactions between the past and the present of all processes are described in terms of linear regressions. This formulation presupposes that the considered model captures the entire interaction structure of the observed processes. If this is not the case, the causal interpretation provided by the regression coefficients may not be justified, and the GC may measure misleading patterns of causality. The classical strictly causal VAR model commonly adopted for GC estimation is limited in this respect, because its coefficients describe the time-lagged effects between the processes but do not account for the instantaneous (i.e. not lagged) effects. As a consequence of this incomplete description, the presence of significant instantaneous effects may lead the standard VAR models to provide wrong interpretations about causality. In fact, recent studies (Hyvarinen et al 2010, Faes 2014) demonstrated the adverse impact of neglecting instantaneous effects on causality analysis based on VAR modeling, providing theoretical arguments showing that the omission of zero-lag correlations from the model can change drastically the values of time-lagged coefficients, and thus of the causality measures based on the model. As a possible solution, these studies proposed the utilization of an extended model accounting for both instantaneous and lagged effects. The extended model was devised combining a classical VAR model and a model of the instantaneous causal interactions among the observed variables commonly known as structural equation model (SEM). Model estimation was performed using classic least-squares regression for the VAR part, and a more sophisticated method called LiNGAM (linear non-Gaussian acyclic model) (Shimizu et al 2006) which exploits non-Gaussianity to address the known identifiability problem of SEM models. The exploitation of combined VAR and SEM models for the study of a generalized form of GC which aggregates both instantaneous and time lagged-effects started with the works 828
L Schiatti et al
Physiol. Meas. 36 (2015) 827
of Hyvarinen et al (2008, 2010), and was employed later on to propose improved measures of causality in the frequency domain (Faes et al 2013a). In the present work we focus on the impact of neglecting instantaneous effects on the computation of the traditional time-domain measure of GC from X to Y, that is defined relating the prediction error variances of two linear regressions, respectively including and excluding the past of X in the prediction of the present of Y. This approach is ubiquitously used in time series analysis without considering the instantaneous effects among the modeled time series (Brovelli et al 2004, Wang et al 2007, Barnett and Seth 2014). Here we show that, since GC reflects the strictly causal VAR representation of the interactions, computing it when instantaneous effects are not trivial may lead to wrong conclusions about the true Granger-causal effects. We propose a measure of extended GC (eGC) that overcomes the limitations of GC by properly accounting for zerolag effects in the linear regressions. We also present a procedure for the practical estimation of eGC which addresses the problem of assigning instantaneous effects to the regression models. This procedure employs a two-step approach based on first estimating the existence of zero-lag correlations between two processes, and then orienting them along one of the two possible causal directions using pairwise measures of non-Gaussianity (Hyvarinen and Smith 2013). After validating the proposed method in a simulation study, we applied it to a practical case of study: the characterization of cardiovascular (CV) and cerebrovascular (CB) regulation mechanisms in subjects prone to postural-related syncope (PRS). PRS is an important and common clinical problem, consisting in a transient loss of consciousness and postural tone, with spontaneous recovery. Its underlying pathophysiological mechanisms are still not fully understood, although it has been associated with an impairment of CV and CB regulation (Toyry et al 1997, Dan et al 2002, Bechir et al 2003, Faes et al 2013a, 2013b). In this study we observe CV and CB variability measuring the fluctuations of heart period (HP), arterial pressure (AP), and cerebral blood flow velocity (FV) during a head-up tilt (HUT) testing protocol that is adopted in clinics to evoke PRS. Then, we compute traditional and novel measures of GC on the three time series measured during specific epochs of the test, in order to determine how the patterns of directional interaction between these variables are affected by the development of PRS. In this physiological application, the extension of the analysis to the inclusion of instantaneous effects is motivated by the fact that zero-lag correlations are expected to be relevant as they are related to fast (i.e. within-beat) mechanisms of CV and CB regulation. In CV regulation, within-beat interactions between systolic AP and HP are compatible with the short latencies, typically lower than 1 s, of the human cardiac baroreflex (Borst and Karemaker 1983); the presence of these interactions has been associated with the fast response of the parasympathetic component of the autonomic nervous system, and has been studied during postural stress and PRS (Westerhof et al 2006, Faes et al 2013b). In CB regulation, studies of dynamic cerebral autoregulation have shown that the response of cerebral blood FV to fast transient changes in AP is completed within a few seconds, and occurs in large part before the next cardiac pulsation (Zhang et al 1998, Panerai et al 1999). The paper is structured as follows: in section 2 the proposed framework for the eGC estimation based on VAR modeling is presented; in section 3 a simulation study to test the effectiveness of eGC measure when instantaneous effects are present is shown; in section 4 the application of the proposed method to a real data set is described; section 5 concludes the paper. We distribute a Matlab toolbox for the computation of GC and eGC; the toolbox is available online at www.lucafaes.net/eGC.html.
829
L Schiatti et al
Physiol. Meas. 36 (2015) 827
2. Methods 2.1. Extended GC
GC is assessed in its traditional formulation based on VAR models (Bressler and Seth 2011). Specifically, given a multivariate stochastic process Y composed of M scalar processes Y1,...,YM, their dynamical interactions are described by the VAR model of order p p
Y (n ) =
∑ A (k ) Y (n − k )
+ U (n ) ,
(1)
k=1
where Y(n) = [Y1(n)⋅⋅⋅YM(n)]T is the vector variable obtained sampling the process Y at the time instant n, A(k) are M × M coefficient matrices in which the element Aji(k) describes the dependence of Yj(n) on Yi(n − k) (i,j = 1,...M; k = 1,...,p), and U = [U1⋅⋅⋅UM]T represents the model residuals. To compute GC from the ith process Yi to the jth process Yj, two linear regressions are considered starting from the VAR representation in (1). First, Yj is described in terms of the past of all processes in Y by means of a so-called unrestricted regression that is obtained considering the part of the joint representation in (1) relevant to the jth process p
Y j (n ) =
∑ Aj (k ) Y (n − k )
+ Uj (n ) ,
(2)
k=1
where Aj(k) is the jth row of A(k). Then, Yj is described in terms of the past of all processes in Y except Yi by means of the so-called restricted regression p
Y j (n ) =
∼
∑ Aj (k ) Y i (n − k )
∼ + Uj (n ) ,
(3)
k=1
where Yi = [Y1⋅⋅⋅Yi–1Yi+1⋅⋅⋅YM]T denotes the original vector process devoid of the ith scalar ∼ ∼ process. Note that the coefficients Aj and the residuals Uj yielded by the restricted regression (3) are in general different than the Aj and the Uj yielded by the unrestricted regression (2). The GC measure is computed by comparing the variance of the residuals resulting from the two regressions: ∼ var Uj (n ) GC . i → j = ln (4) var Uj (n ) The GC measure (4) is always positive and proportional to the extent to which the knowledge of the past of Yi helps improving the prediction of Yj(n) above and beyond the improvement achieved using the past of Yi; the measure becomes null if the past of Yi is ineffective to the prediction improvement (in this case, the coefficients Aji(k) are zero for all lags k and the restricted and unrestricted regressions (2) and (3) coincide). An important property of the GC is that it is sensitive only to the causal effects from Yi to Yj which are not mediated by the interaction of Yi and Yj with any other process Yk (k ≠ i, k ≠ j). Therefore, the GC is able to elicit direct dependencies between two processes and to discard the confounding effects of third processes belonging to the set of processes used for VAR modeling. The model of equation (1) is a strictly causal model, in the sense that it describes only the time-lagged interactions between the processes. Therefore, GC is computed without considering the instantaneous (i.e. non-delayed) effects among the observed time series. However, the possible presence of zero-lag effects is known to have an impact also on the time-lagged effects (Hyvarinen et al 2010, Faes 2014), and thus may affect the reliability of the observed 830
L Schiatti et al
Physiol. Meas. 36 (2015) 827
GC patterns. To overcome this problem, we propose the utilization of an extended version of GC, derived from an extended VAR model combining both instantaneous and lagged effects: p
Y (n ) =
∑ B (k ) Y (n − k )
+ W (n ) ,
(5)
k=0
where B(k) are M × M coefficient matrices in which the element Bji(k) describes the dependence of Yj(n) on Yi(n − k) (i,j = 1,...M; k = 0,...,p), including also coefficients weighing the present samples of Y, and W = [W1⋅⋅⋅WM]T represents the model residuals. To compute eGC from Yi to Yj, unrestricted and restricted linear regressions incorporating zero-lag causal effects are computed starting from the extended VAR representation in (5), analogously to what seen for GC: p
Y j (n ) =
∑ Bj (k ) Y (n − k )
+ Wj (n ) ,
(6)
+ W͠ j (n ) ,
(7)
k=0 p
Y j (n ) =
∼
∑ Bj (k ) Y i (n − k ) k=0
∼ where Bj(k) and Bj (k ) are row coefficients weighting respectively the present and past samples of all processes in Y for the unrestricted regression (6), and the present and past samples of all processes in Y except Yi for the restricted regression (7). The eGC measure is computed by comparing the variance of the residuals of the two regressions: eGC i → j = ln
var W͠ j (n ) var Wj (n )
.
(8)
2.2. Estimation
In practical analysis, GC and eGC are estimated between any pair of time series taken from an observed multivariate dataset by identifying the unrestricted and restricted regression models (respectively, equations (2) and (6) and equations (3) and (7)) and then combining the variance of the estimated model residuals as in equations (4) and (8). In this study, regressions were performed by the standard least squares approach, and the number of lags p was set identical for all time series by applying the Bayesian Information Criterion to the overall VAR representation of equation (1) (Faes et al 2012). Moreover, the statistical significance of each computed GC and eGC measure was assessed using the Fisher F-test. While for the computation of the traditional GC the estimation is straightforward as just described above, for the eGC analysis the procedure is complicated by the need of identifying the structure of the zero-lag effects among the time series. Such a structure cannot be imposed arbitrarily, because it has to be consistent with the joint representation of the observed multivariate time series (Shimizu et al 2006). The zero-lag effects are collected in the matrix B(0) of the extended VAR representation in (5): allowing an effect at lag zero between two processes Yi and Yj corresponds to set either Bij(0) ≠ 0 or Bji(0) ≠ 0; the coefficients cannot be both nonzero because the presence of a pairwise loop would make the extended model unidentifiable (Shimizu et al 2006, Hyvarinen et al 2010, Faes et al 2013a). Thus, the estimation of eGC requires to identify the presence and causal direction of the zero-lag effects between the observed processes to be incorporated into the regression models (6) and (7). In the following, we propose a two-step procedure to accomplish this task, based first on estimating the 831
L Schiatti et al
Physiol. Meas. 36 (2015) 827
existence of zero-lag correlations in an undirected sense, and then on finding their directions using pairwise measures of non-Gaussianity (Hyvarinen and Smith 2013). As a first step, we determine whether an effect at lag zero is present between any pair of processes Yi and Yj, looking at the correlation between the variables Yi(n) and Yj(n). The simplest way to assess correlations in a multivariate process is to compute its partial correlation matrix, which is defined as the inverse of the covariance matrix of the process (Marrelec et al 2006). Accordingly, in our preliminary study (Schiatti et al 2014) we computed the partial correlation of Y, [covY(n)]–1. However, the partial correlation reflects correctly the true correlations in the data only in the absence of any temporal structure, while it may lead to detect spurious instantaneous correlations when lagged correlations are also present. This occurs due to possible common driver effects from some past lag component, say Yk(n − l), which drives both Yi(n) and Yj(n) inducing a zero-lag correlation that is not conditioned out by the partial correlation. On the other hand, assuming that the VAR model (1) properly captures all time-lagged effects into its coefficients, the residuals U do not have any temporal structure, so that finding nonzero partial correlation between Ui(n) and Uj(n) means that a ‘true’ direct effect exists between Yi(n) and Yj(n). Therefore in this study we assessed zero-lag correlations in our multivariate process computing the partial correlation of the VAR model residuals, [covU(n)]−1. In the practical computation, significant nonzero values of the estimated partial correlations were determined using bootstrap resampling (Zoubir and Boashash 1998). Specifically, we drew 100 bootstrap samples of the same length from the original estimates of U(n) through sampling with replacement and computed the partial correlation for each sample. We then exploited the resulting distribution of values to get a 95% confidence interval for [covU(n)]−1, and labeled as statistically significant all elements of the original partial correlation matrix for which the confidence interval did not encompass the zero level. In the second step, the direction of the zero-lag interaction between pairs of series associated with residuals showing a significant partial correlation was assessed employing the pairwise likelihood ratios for non-Gaussian data recently proposed by Hyvärinen and Smith (2013). This approach is essentially based on estimating the likelihood of the models corresponding to the two directions of causal influence between a pair of variables, and on deciding in favor of the model with higher likelihood after the computation of a likelihood ratio. The two models, which consider the standardized (i.e. having zero mean and unit variance) variables x and y, are defined by the simple regressions y = ρx + d (model x → y) and x = ρy + e (model y → x), where d and e are the residuals, assumed to be independent of the regressors x and y, and ρ is the correlation coefficient between x and y, and is therefore the same in the two models. Considering N realizations of x and y, the ratio between the likelihoods of the two models, L(x → y) and L(y → x), is first computed as R =
L (x → y ) 1 ; log L (y → x ) N
(9)
then, the causal direction is chosen as x → y if R is positive, and as y → x if R is negative. In (Hyvarinen and Smith 2013), exploiting the information theoretic interpretation of the likelihood ratio and the connection between information measures and non-Gaussianity, it is shown that choosing the model with higher likelihood corresponds to choosing the direction which leads to maximum non-Gaussianity. This is in line with the knowledge that the direction of SEM models describing the instantaneous interactions among processes can be inferred only if the observed variables have a non-Gaussian distribution (Shimizu et al 2006, Hyvarinen et al 2010, Faes et al 2013a). For the practical computation of the likelihood ratio, we adopted the first-order approximation of the likelihood of the two models, describing the 832
L Schiatti et al
Physiol. Meas. 36 (2015) 827
log-probability density functions with the logistic density (Hyvarinen and Smith 2013). This leads to estimating the likelihood ratio as R ̂ = ρ ̂ x tanh(y )− y tanh(x ) ,
(10)
where ρ ̂ is the estimated correlation coefficient and denotes sample average. In our application, the likelihood ratio was computed as in (10) taking as variables x and y the standardized version of any pair of estimated VAR residuals Ui(n) and Uj(n) for which a statistically significant partial correlation was found according to bootstrap resampling. 3. Simulation study 3.1. Theoretical formulation
In order to describe the features of GC and eGC, and to illustrate their relationship, we consider the VAR process of order p = 2, composed by M = 3 scalar processes, generated by the equations:
⎧Y1(n )= δY2(n ) + b11(1)Y1(n − 1) + b11(2)Y1(n − 2) + W1(n ) ⎪ ⎨Y2(n )= b 22(1) Y2(n − 1) + b 22(2) Y2(n − 2) + b 21(2)Y1(n − 2) + b 23(2)Y3(n − 2) + W2(n ) , (11) ⎪Y (n )= δY (n ) + b (1)Y (n − 1) + W (n ) ⎩ 3 1 31 1 3
where the Wi are white and uncorrelated model residuals. The coefficients appearing on the right-hand side of the equations can be collected in the matrices B(k), k = 0,1,2. The diagonal elements B(k) are set in order to generate complex conjugate poles for the processes Y1 and Y2 (b11(1) = 2r1cos(2πf1); b11 (2) = −r12; b22(1) = 2r2cos(2πf2); b 22 (2) = −r22); we set modulus r1 = 0.9 and r2 = 0.8, and phase ±3/5π and ±π/5 determining autonomous oscillations at f1 = 0.3 Hz and f2 = 0.1 Hz, respectively for Y1 and Y2. Then, the off-diagonal elements of B(k) determine direct causal effects from one process to another. The causal effects at lag zero are set according to the parameter δ; here we consider the conditions δ = 0 and δ = 0.8, entailing absence and presence of instantaneous effects, respectively. The process generated by (11) is an extended VAR process, which can be described in a compact form as in equation (5). The corresponding strictly causal VAR representation can be obtained analytically from the extended representation using the relations (11): U (n ) = [I − B(0)]−1 W (n ),
(12)
A(k ) = [I − B(0) ]−1 B(k ) for each k ≥ 1.
(13)
These relations allow to compute the model coefficients and residuals of the strictly causal VAR model equivalent to the extended model defined by (11). Doing so, we depicted the patterns of causality resulting from the two representations in the two conditions δ = 0 and δ = 0.8, obtaining the connectivity diagrams of figure 1. When δ = 0 all the imposed effects from one process to another are lagged, and result from the nonzero model coefficients: Y1 → Y2 (b21(2) = 0.5), Y1 → Y3 (b31(1) = 0.5) and Y3 → Y2 (b23(2) = 0.5). In this condition with absence of instantaneous causality the extended and strictly causal VAR representations coincide (B(0) = 0 in (12) and (13) makes U(n) = W(n) and A(k) = B(k)), so that both models provide a correct interpretation of the causality patterns (figure 1(a)). When δ = 0.8, two instantaneous effects over the directions Y2 → Y1 and Y1 → Y3 are imposed in addition to the lagged causal effects described so far. In this case the extended VAR representation correctly reflects relations imposed by equation (11), as depicted in figure 1(b) (left). On the contrary, due to the 833
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Figure 1. Extended VAR representation (left) and strictly causal VAR representation (right) of the VAR processes generated by the equations in (11), in absence (a) and presence (b) of instantaneous effects. Arrows represent causal effects, with the numbers indicating the lag at which effects are detected.
presence of instantaneous causality the causal interpretation based on the strictly causal model is now different, as it reflects the strictly causal model coefficients computed using (13). In particular, we see that the strictly causal representation identifies more lagged relations than the imposed ones, interpreting as lagged the instantaneous effect Y2 → Y1 and detecting the spurious lagged causality effects Y3 → Y1 and Y2 → Y3, as well as an unexpected effect from Y3(n − 2) to Y3(n). This example shows that the traditional strictly causal representation may provide misleading patterns of causality when the instantaneous effects among the observed processes are not negligible. 3.2. Implementation on simulations
The feasibility of the estimation approaches illustrated in section 2.2 for the computation of GC and eGC measures was tested on simulations of the extended VAR process illustrated in section 3.1. Specifically, realizations of (11) lasting N = 300 points were generated starting from realizations of the extended residuals Wi, and the comparative analysis was performed under three different scenarios: (a) absence of instantaneous effects (δ = 0 in (11)); (b) presence of instantaneous effects (δ = 0.8) with non-Gaussian extended residuals; and (c) presence of instantaneous effects (δ = 0.8) with Gaussian extended residuals. The extended residuals, Wm(n), m = 1,…,M, n = 1,…,N, were generated as independent Gaussian white noises for scenarios (a) and (c), and as independent non-Gaussian white noises for scenario (b); in this last case, non-Gaussianity was achieved first generating Zm(n) as independent Gaussian white noises and then applying the nonlinear transformation Wm(n) = sign(Zm(n))|Zm(n)|q, with exponent q chosen in the range [0.5, 0.8] or [1.2, 2.0] to yield respectively a sub-Gaussian or super-Gaussian distribution for Wi (Faes et al 2013a). GC and eGC were computed from 100 realizations of the simulation generated as described above for the three scenarios. Moreover, the theoretical values of GC and eGC were also computed exploiting the procedure described in (Faes et al 2015, Barnett and Seth 2014). This procedure allows to express the variance of the residuals of both unrestricted and restricted models, and thus to express the GC and eGC, in terms of the VAR parameters. Here, theoretical 834
L Schiatti et al
Physiol. Meas. 36 (2015) 827
values of GC and eGC were used to assess the reliability of the measures estimated from the process realizations. Figure 2 reports the distributions of GC and eGC estimated over 100 realizations of the simulation under the three scenarios described above. Moreover, table 1 reports the number of realizations for which instantaneous effects were detected along any specified causal direction. In the absence of instantaneous effects, the distributions of the estimated GC and eGC overlap almost perfectly, and are centered in the theoretical values (figure 2(a)). Moreover, both estimates were always detected as statistically significant, according to the F-test, when computed along the coupled directions Y1 → Y2, Y1 → Y3 and Y3 → Y2, and were almost never detected as significant when computed along the uncoupled directions. This demonstrates that both approaches are able to retrieve the patterns of GC when the investigated processes do not interact at lag zero. The equivalence between GC and eGC in this case is confirmed by the observation that the bootstrap procedure performed to test for instantaneous causality detected the presence of significant partial correlations among the VAR residuals in a very limited number of realizations (table 1(a)). When significant instantaneous effects were simulated imposing δ = 0.8 in (11), the estimates of GC and eGC were still centered on their theoretical values (figure 2(b)), documenting the appropriateness of the estimation procedures. However, in this case the patterns of GC and eGC were substantially different: eGC was markedly higher than GC along the directions Y2 → Y1 and Y1 → Y3, and markedly lower and non-significant along the directions Y3 → Y1 and Y2 → Y3. Comparing figure 2(b) with figure 1(b), we see that the two measures reflect the extended and strictly causal representations of the simulated processes: the eGC is significant whenever an instantaneous or time lagged effect is set from one process to another, while the GC follows the misleading patterns of lagged causality provided by the strictly causal representation. This confirms the unsuitability of GC to measure Granger-causal effects in processes exhibiting zero-lag correlations. On the contrary, the estimation procedure based on the extended VAR representation correctly detects the presence and direction of zero-lag interactions among the processes (table 1(b)) and provides a meaningful causal interpretation based on the eGC (figure 2(b)). Figure 2(c) illustrates what happens when the requirement of non-Gaussianity of model residuals is not met. In this case, instantaneous effects can be detected if not present, or be assigned to the wrong causal direction (table 1(c)). As a consequence the distributions of the estimated eGC became wider and deviated from the theoretical values (figure 2(c)). This documents the importance of satisfying the requirement of non-Gaussianity to guarantee a proper identification of instantaneous effects, and consequently a reliable estimation of Grangercausal effects. 4. Real data application This section describes the application of the proposed framework based on VAR modeling to estimate patterns of GC and eGC in physiological time series measured in the context of CV and CB regulation in subjects with PRS. 4.1. Experimental protocol
Ten young subjects (3 males, 22.4 ± 9.5 years) with history of recurrent unexplained syncope were studied during a HUT test protocol. The acquired signals were the electrocardiogram (ECG), the finger photoplethysmographic AP, and the cerebral blood flow based on 835
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Figure 2. Box-plot of GC and eGC values computed over 100 realizations of equation (11) generated according to the three scenarios: (a) δ = 0, Gaussian Wi; (b) δ = 0.8, non-Gaussian Wi; (c) δ = 0.8, Gaussian Wi. The symbols ✠ indicate the theoretical value of GC and eGC. The number of GC or eGC values detected as statistically significant (F-test with 0.01 significance) is reported below each distribution.
transcranial Doppler ultrasonography, measured with subjects lying in the resting supine position and standing in the 60° upright position reached after passive HUT. Postural stress in the upright position was prolonged until the occurrence of PRS. All subjects exhibited spontaneous recovery with return to the supine position. The analyzed time series were the beat-to-beat variability of the HP, the mean AP and the mean cerebral blood FV, measured respectively as the temporal distance between two consecutive R peaks of the ECG, and as the integral of the AP and cerebral flow signals taken within each detected diastolic pulse interval normalized to interval duration (figure 3). Three stationary epochs of N = 300 beats were considered for the analysis, one in the supine position (SU), and two in the early (ET) and late (LT) phases of passive standing in the upright position (respectively at ~2 min after tilting and immediately before the drop in AP indicating presyncope). Details about experimental protocol can be found in (Faes et al 2013c). 836
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Table 1. Percentage of realizations of equation (11) for which instantaneous correlations
were detected along any given causal direction, looking at the statistical significance of the partial correlation of the VAR residuals, under the three scenarios: (a) δ = 0, Gaussian Wi; (b) δ = 0.8, non-Gaussian Wi; (c) δ = 0.8, Gaussian Wi. Zero-lag effects
(a) (b) (c)
Y2 → Y1
Y1 → Y2
Y3 → Y1
Y1 → Y3
Y3 → Y2
Y2 → Y3
6 100 52
1 0 48
1 0 46
2 100 54
3 0 4
3 8 5
Figure 3. HP, AP and FV variability series measure; (a) computation of HP from ECG signal; (b) computation of AP from finger photoplethysmographic AP signal; (c) computation of FV from transcranial Doppler ultrasonography signal.
4.2. Data analysis
GC and eGC were computed between each pair of the three considered time series in the SU, ET and LT windows. The model order was assessed using the Bayesian information criterion (Faes et al 2012), minimizing the corresponding figure of merit in the range 1–20. The identification of zero-lag correlations was performed using 100 bootstrap samples and a 95th percentile confidence interval. The F-test for assessing nonzero values of GC and eGC was run with 1% significance. The statistical significance of the differences in the values of GC or eGC across the three conditions (SU, ET, LT) was assessed by the one-way ANOVA with 5% significance, followed by post-hoc comparison based on the Student t-test for paired data with Bonferroni correction. Moreover, the statistical significance of the difference between the distributions of GC and eGC computed for a specific condition was assessed by the Student t-test for paired data. In addition, since the identification of instantaneous effects relies on the hypothesis of nonGaussianity for the model residuals, such an hypothesis was tested by means of a method based on mutual information (MI) (Hlinka et al 2011). This method exploits the knowledge that the MI I(x,y) between two standardized variables x and y has a lower bound such that I(x,y) ≥ IGauss(r) 837
L Schiatti et al
Physiol. Meas. 36 (2015) 827
= −0.5ln(1−r2), where r is the linear correlation between x and y; the equality holds exactly for a bivariate Gaussian distribution, for which the MI IGauss is a function of the correlation between the two variables. This allows to assess the deviation from Gaussianity of the joint distribution between x and y by comparing the MI I(x,y) with the MI between surrogates of x and y which have the same linear correlation of the original variables but are enforced to have bivariate Gaussian distribution. In our study, the role of x and y was taken by any pair of estimated VAR residuals Ui(n) and Uj(n), and surrogates were constructed as multivariate Fourier transform (FT) surrogates (Prichard and Theiler 1994). These surrogates were obtained computing the FT of the original series, keeping the magnitudes of the Fourier coefficients but adding the same random number to their phase, and finally returning into the time domain through inverse FT computation. We generated 100 surrogate pairs for each original pair of residuals, and rejected the null hypothesis of two linear Gaussian processes when the original MI was larger than the 95th percentile of the distribution of MI computed over the surrogate pairs. The MI was computed using the nearest neighbor estimator proposed by Kraskov et al (2004). 4.3. Results and discussion
The order of the regression models employed for the computation of GC and eGC, selected with the Bayesian Information Criterion, was comprised for all subjects between 2 and 5 and did not vary significantly across the three considered epochs of HUT (supine: 3.6 ± 1.1; early tilt: 3.1 ± 0,6; late tilt: 3.4 ± 0.5 mean ± SD across subjects). The analysis of GC and eGC, performed in terms of both magnitude of the measure and number of subjects with significant causality, led to the results depicted in figure 4. This analysis indicated the existence of a bidirectional causality between HP and AP, with prevalence of HP → AP causality, and of mostly unidirectional causality along the directions HP → FV and AP → FV. Moreover we found significantly higher causality AP → HP during ET compared with SU, and higher causality AP → FV during LT compared with SU. These results are in line with expected mechanisms of CV and CB regulation. Indeed, the bidirectional HP–AP interaction reflects the closed loop regulation of the CV system, resulting from the coexistence between the feedback baroreflex mechanism acting from AP to HP and the feedforward mechanical effects of HP on AP (Cohen and Taylor 2002). The strengthening of the causal coupling from AP to HP with the tilt transition documents the well known activation of the baroreflex (Porta et al 2011, Faes et al 2013b), developed as a reaction to the change of posture finalized to maintain blood pressure to a constant level; on the opposite we observed a decrease in the baroreflex coupling AP → HP moving from ET to LT, though not statistically significant, that might reflect an impairment of this mechanism in proximity of syncope (Faes et al 2006). As to the analysis of CB regulation, the significant causal effect of the two CV variables on FV are expected since it is known that the cerebral blood flow can vary rapidly in response to systemic AP variations (Zhang et al 2000). On the other hand, it is also known that the CB bed works for maintaining a relatively constant blood flow in the face of the perturbing action of AP, according to a homeostatic mechanism known as cerebral autoregulation (Paulson et al 1990). According to this mechanism, increases in the coupling between AP and FV have been often associated with a loss of cerebral autoregulation (Panerai et al 1998, Mitsis et al 2004), and in this view the higher GC from AP to FV which we observed upon presyncope may be indicative of impaired CB control. The results described above were confirmed using eGC, which indeed increased from SU to ET in the direction AP → HP (though without reaching statistical significance, p = 0.114), and increased significantly from SU to LT in the direction AP → FV. Furthermore, eGC analysis showed modifications of the CB causality patterns associated with presyncope that were 838
L Schiatti et al
Physiol. Meas. 36 (2015) 827
Figure 4. Box-plot of GC and eGC values between HP, mean AP and mean FV series computed during SU, ET and LT epochs of HUT. The number of GC or eGC values detected as statistically significant (F-test with 0.01 significance) is reported below each distribution. Significant difference between pairs of distributions: *p
