A new Granger causality measure for eliminating the confounding influence of latent common inputs

TL;DR

This paper introduces a new Granger causality measure that effectively removes the confounding effects of latent common inputs, improving the accuracy of detecting true directed interactions in time series data.

Contribution

A novel Granger causality measure inspired by partial Granger causality that is robust against latent confounders, with a practical inference procedure demonstrated through numerical experiments.

Findings

The test statistics follow F-distributions under no interaction.

The proposed method eliminates false positives caused by latent confounders.

Normal Granger causality detects spurious interactions due to confounders.

Abstract

In this paper, we propose a new Granger causality measure which is robust against the confounding influence of latent common inputs. This measure is inspired by partial Granger causality in the literature, and its variant. Using numerical experiments we first show that the test statistics for detecting directed interactions between time series approximately obey the $F$-distributions when there are no interactions. Then, we propose a practical procedure for inferring directed interactions, which is based on the idea of multiple statistical test in situations where the confounding influence of latent common inputs may exist. The results of numerical experiments demonstrate that the proposed method successfully eliminates the influence of latent common inputs while the normal Granger causality method detects spurious interactions due to the influence of the confounder.

Tables

Table 1. Table 1: The values of the model parameters for each simulated model

model	$a_1$	$a_2$	$b_1$	$b_2$	$c_1$	$c_2$	${\sigma }_x$	${\sigma }_y$	$\rho$
$\left(a\right)$	0.9	$-0.5$	0.5	$-0.2$	0	0	1	$\sqrt{0.7}$	0
$\left(b\right)$	0.9	$-0.5$	0.5	$-0.2$	0	0	1	$\sqrt{0.7}$	0.4
$\left(c\right)$	0.9	$-0.5$	0.5	$-0.2$	0.16	$-0.2$	1	$\sqrt{0.7}$	0
$\left(d\right)$	0.9	$-0.5$	0.5	$-0.2$	0.16	$-0.2$	1	$\sqrt{0.7}$	0.4

Equations

$$X_t=\sum _{i=1}^{\infty }{a}_{1i}{X}_{t-i}+\sum _{i=1}^{\infty }{c}_{1i}{Z}_{t-i}+u_{1t},$$

$$Z_t=\sum _{i=1}^{\infty }{b}_{1i}{Z}_{t-i}+\sum _{i=1}^{\infty }{d}_{1i}{X}_{t-i}+u_{2t},$$

$$S=$$

$$=$$

$$X_t=\sum _{i=1}^{\infty }{a}_{2i}{X}_{t-i}+\sum _{i=1}^{\infty }{b}_{2i}{Y}_{t-i}+\sum _{i=1}^{\infty }{c}_{2i}{Z}_{t-i}+u_{3t},$$

$$Y_t=\sum _{i=1}^{\infty }{d}_{2i}{X}_{t-i}+\sum _{i=1}^{\infty }{e}_{2i}{Y}_{t-i}+\sum _{i=1}^{\infty }{f}_{2i}{Z}_{t-i}+u_{4t},$$

$$Z_t=\sum _{i=1}^{\infty }{g}_{2i}{X}_{t-i}+\sum _{i=1}^{\infty }{h}_{2i}{Y}_{t-i}+\sum _{i=1}^{\infty }{k}_{2i}{Z}_{t-i}+u_{5t},$$

$$\Sigma =$$

$$=$$

$$R_{XX|Z}^{(1)}\equiv$$

$$R_{XX|Z}^{(2)}\equiv$$

$$\displaystyle F_{1}=\ln\Biggl{(}\frac{R_{XX|Z}^{(1)}}{R_{XX|Z}^{(2)}}\Biggr{)}.$$

$$x_t=$$

$$y_t=$$

$$\Sigma =$$

$$=$$

$$\Sigma =\left(\begin{array}{cc}{\tau }^2+{\eta }^2{\sigma }_y^2& \eta {\sigma }_y^2\\ \eta {\sigma }_y^2& {\sigma }_y^2\end{array}\right),$$

$${\tau }^2={\Sigma }_{xx}-{\Sigma }_{xy}{\Sigma }_{yy}^{-1}{\Sigma }_{yx}$$

$${\epsilon }_t=\eta {\xi }_{t-1}+{\omega }_t,$$

$${\omega }_t\sim \mathcal{N}(0,{\tau }^2),$$

$$\mathrm{cov}({\omega }_t,{\xi }_{t-1})=\mathrm{cov}({\omega }_t,{\epsilon }_t)=0.$$

$$x_t=$$

$$y_t=$$

$$\{\begin{array}{ll}{H}_0:& c=0,\\ H_1:& c\ne 0.\end{array}$$

$$p(x_t|x_{t-1},y_{t-1},y_{t-2})=\mathcal{N}({\mu }_t^{(1)},{\Sigma }_{x|y}),$$

$${\mu }_t^{(1)}=$$

$$=$$

$${\Sigma }_{x|y}=$$

$$=$$

$$L(\mathbf{\beta },b|\mathbf{x},\mathbf{y})=p(x_T,x_{T-1},\dots ,x_3,y_{T-1},\dots ,y_1|x_2,y_1,\mathbf{\beta },b).$$

$$L(\mathbf{\beta },b|\mathbf{x},\mathbf{y})=$$

$$=$$

$$⋮$$

$$=$$

$$\equiv$$

$$\log L_{x|y}(\mathbf{\beta },b|\mathbf{x},\mathbf{y})=$$

$$=$$

$$\equiv$$

$$\log L_y(b|\mathbf{y})=$$

$$D_1\equiv$$

Full text

I Introduction

The Granger causality test is a method for inferring causal interactions between time series by a statistical test. This method is based on the idea by Wiener that if past information of one time series improves the prediction of another time series, there is a causal influence Wiener (1956). Granger formalized this notion of causality in the framework on autoregressive (AR) modeling of multivariate time series in the econometrics literature Granger (1969). Strictly speaking, Granger causality does not mean causality in the true sense, although it should be considered as a necessary condition for the true causality. However, the method has gained popularity in a wide range of fields, such as econometrics and biology due to its conceptual simplicity and ease of implementation Granger (1980); Baccala et al. (1998); Gourevitch et al. (2006). For example, in the field of biology, the Granger causality test was applied to electroencephalogram Barrett et al. (2015), functional magnetic resonance imaging of Blood Oxygen Level Dependent signals Ding et al. (2006); Chen et al. (2006); Bressler and Seth (2011); Seth et al. (2015), multi-electrode arrays of local field potentials and Calcium imaging of neuronal activity De Vico Fallani et al. (2015). Also, more recently, Granger causality has extended to point process time series and applied to spike activity of neurons Stevenson et al. (2009); Kim et al. (2011); Gerhard et al. (2013).

Nowadays, application of the Granger causality test is more common. However, when we apply the Granger causality test to real data, we must pay attention to confounding environmental influences. In fact, in many experimental settings, we are only able to record a subset of all related variables in a system. In such a case, applications of Granger causality may detect spurious directed interactions due to the influence of latent common inputs and unobserved variables. For example, one time series may falsely appear to cause another if they are both influenced by a third time series or a common input but with a delay. In addition, data analysis in practice involves the step of model selection, in which a relevant set of variables is selected for analysis Roebroeck et al. (2011). Then, this step is likely to exclude some relevant variables, which can lead to the detection of apparent causal interactions that are actually spurious Pearl (2000). Hence, controlling for latent common inputs and latent variables is a critical issue when applying Granger causality to experimental data.

Partial Granger causality was developed as a method for eliminating the confounding influence of latent common inputs and unobserved variables Guo et al. (2008). According to what the authors in the literature say, the method was inspired by the partial correlation in statistics. They insisted that partial Granger causality does estimate directed interactions more robustly than conditional Granger causality against the influence of latent common inputs and latent variables. In the recent literature we can find some reviews and applications of this type of Granger causality Smith et al. (2011). Although the idea of partial Granger causality is fascinating and the method is applied to various systems, the justification of the method is arguable Roelstraete and Rosseel (2012). In the original paper, it was insisted that the sampling distribution of the test statistics $F_{1}$ cannot be determined analytically, and one has to resort to the bootstrap technique. Furthermore, the test statistics can take a negative value, which was reasoned as being subtracted by an additional term. By contrast, in the subsequent study in which partial Granger causality was extended to multivariate system, it was argued that the test statistics cannot take a negative value Barrett et al. (2010). In another literature, it was argued that the negative value in the test statistics of partial Granger causality is a serious flaw, and this undermines the credibility of the obtained results and, thus, the validity of the approach Roelstraete and Rosseel (2012). Therefore, it remains unclear that at which circumstance partial Granger causality is available and not even in numerical experiments.

In this paper, modifying the idea of partial Granger causality, we propose a new Granger causality measure which is robust against the confounding influence of latent common inputs. We first point out that the original partial Granger causality is attempting to eliminate the influence of the confounder between the target variable and conditioning variable, rather than the target variable and source variable. Then, we introduce a Granger causality measure by eliminating the influence of latent common inputs between the target variable and source variable. We also discuss the way of parameter estimation for the model. In our parameter estimation procedure, the test statistics cannot take a negative value. Of the influence of latent confounder, this paper only considers the latent common inputs, and the latent variables are not considered which has a worse effect on the estimation of directed interactions Roelstraete and Rosseel (2012). The reason for this is that the influence of latent variables have to be dealt with in the framework of the Moving Average (MA) model, rather than the AR model.

This paper is organized as follows. In Sec. II, we review partial Granger causality and introduce our Granger causality measure to illustrate and highlight the difference between partial Granger causality and proposed Granger causality introduced in this paper. Furthermore, we discuss the distribution that the test statistics obeys by using the likelihood function of the observed time series. We also discuss the way of parameter estimation for the model and the reason why the test statistics can take a negative value in the literature. In Sec. III, we show the validity of our measure using the data generated by the models, which are variants of the models extensively studied in the literature for numerical experiments. Furthermore, we propose a practical procedure for estimating directed interactions where latent common inputs may exist. Sec. IV is devoted to conclusions.

II Proposed measure for Granger causality

In this section, we review partial Granger causality to illustrate and highlight the difference between partial Granger causality and proposed Granger causality introduced in this paper, and point out the issue of partial Granger causality. Then, we introduce our new Granger causality measure.

II.1 Partial Granger causality

Seth et al., proposed the following idea of Granger causality which eliminates the influence of latent confounder. In their notation, they considered the following AR model for the null hypothesis:

[TABLE]

where $u_{1t}$, $u_{2t}$ are noise terms which have zero mean and involve the influence of latent common inputs and latent variables. They incorporated the influence into the model explicitly by the following covariance matrix:

[TABLE]

In the same way, the model with directed interactions for the alternative hypothesis is expressed as follows:

[TABLE]

where each noise has a following covariance matrix:

[TABLE]

Then, for each model, they introduce the variance of noise for the target variable $X_{t}$ by eliminating the influence of the noise for $Z_{t}$ as

[TABLE]

Since $R_{XX|Z}^{(1)}$, $R_{XX|Z}^{(2)}$ are variance of noise for target variable with the influence of latent common inputs and latent variables eliminated, they insisted that directed interactions under the existence of latent common inputs and latent variables can be estimated by the analysis of variance of $R_{XX|Z}^{(1)}$, $R_{XX|Z}^{(2)}$. The test statistic is given by

[TABLE]

It is because the test statistic $F_{1}$ quantifies the decrease of variance by including interaction terms with the influence of confounder eliminated. They called the above statistical test as partial Granger causality since it was inspired by the partial correlation in statistics.

II.2 Proposed statistical measure

In our opinion, partial Granger causality introduced in the previous subsection has a conceptual issue. The variance, Eqs. (8) and (9), eliminates the noise correlation between the target variable $X_{t}$ and the conditioning variable $Z_{t}$. That is, they eliminate the confounding environmental influence between the target variable and the conditioning variable. However, when we are going to estimate directed interactions from the source variable $Y_{t}$ to the target variable $X_{t}$, the influence that should be removed is that between the source variable and the target variable, which can cause a spurious interaction from $Y_{t}$ to $X_{t}$. Therefore, we define here the variance by eliminating the influence of latent common inputs between the source variable and the target variable in a simplest example of the AR model. Then, using this variance, we propose a new Granger causality test for directed interactions under the existence of latent common inputs.

The model we suppose consists of two time series of the following AR models $x_{t}$, $y_{t}$, however, the noise structure is different from the usual AR model:

[TABLE]

Here, the noise term $\varepsilon_{t}$ and $\xi_{t}$ are zero mean Gaussian noises, although they have a different-time correlation due to the latent common inputs with delay. That is, the covariance matrix of the noise is expressed as follows,

[TABLE]

where $\rho$ is a correlation coefficient between $\varepsilon_{t}$ and $\xi_{t-1}$, which means that the latent common input reaches $x_{t}$ with delay compared to $y_{t}$. These noises do not have autocorrelation, $\mathrm{cov}(\varepsilon_{t},\varepsilon_{t-i})=\mathrm{cov}(\xi_{t},\xi_{t-i})=0,\hskip 5.69046pti\in\mathbb{N}$. In this paper, we express the influence of latent variables by an auto-correlated noise. Our interest is the estimation of interactions from $y_{t}$ to $x_{t}$, that is, whether the model parameter $c$ exists or not. For simplicity, we assume that the interaction of opposite direction, from $x_{t}$ to $y_{t}$, does not exist, that is, the system of our interest comprises asymmetric network structure.

Here, we note the following fact. In the original paper of partial Granger causality, two types of noise term which are generated by the confounding environmental influence are considered. The first is generated by latent common inputs, which has a cross correlation between each time series but does not have autocorrelation. In addition, only the equal-time cross-correlation was considered. However, the equal-time correlated noise generated by latent common inputs does not produce a spurious directed interaction, this fact will be numerically confirmed later. Latent common inputs with delays could generate a spurious directed interaction when one uses the normal Granger causality test. The influence of latent common inputs can be removed by the method proposed in this paper. The second is the noise generated by the latent variables. The influence of latent variables is expressed by an auto-correlated noise, $\mathrm{cov}(\varepsilon_{t},\varepsilon_{t-i})\neq 0,\hskip 5.69046pti\in\mathbb{N}$. However, strictly speaking, an auto-correlated noise is beyond the scope of the AR model. Therefore, in this paper, we do not suppose the existence of latent variables. The influence of latent variables requires Granger causality constructed by the MA model. The study in this direction has been worked on in the framework of the state space model Barnett and Seth (2015).

To construct our new measure for Granger causality, it is useful to reparametrize the noise parameter $(\sigma_{x},\sigma_{y},\rho)$ to $(\tau,\sigma_{y},\eta)$ as

[TABLE]

where

[TABLE]

is the variance of the noise with the influence of latent common inputs eliminated similar to partial Granger causality. This reparameterization is equivalent to following linear relationship between $\varepsilon_{t}$ and $\xi_{t}$:

[TABLE]

where the noise $\omega_{t}$ has following property:

[TABLE]

[TABLE]

Using the above reparameterization, the model, Eqs. (11) and (12), can be re-expressed as

[TABLE]

Since the conditional variance, $\tau^{2}=\Sigma_{xx}-\Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}$, is a variance of noise for $x_{t}$ from which the correlated component with the noise for $y_{t}$ was eliminated, we propose a statistical measure for Granger causality in detecting directed interactions as whether the variance of $\omega_{t}$ decreases significantly by including the interaction parameter $c$. That is, we set the models of the null hypothesis and the alternative hypothesis as follows:

[TABLE]

II.3 Interpretation by the likelihood function

In this subsection, we re-express the formulation explained in the previous subsection in terms of the likelihood function in order to consider the distribution that our test statistic obeys. First, we note that the variance of the noise, Eq. (15), is formally equivalent to that of the conditional Gaussian distribution. In fact, the probability of the AR model, Eq. (19), can be expressed as a conditional Gaussian distribution:

[TABLE]

where

[TABLE]

Therefore, extracting the correlated component of noise due to latent common inputs is equivalent to conditioning on $y_{t}$.

Now, we consider the situation where time series $x_{t}$ and $y_{t}$ are observed, $\mathbf{x}=(x_{T},x_{T-1},\dots,x_{2})$, $\mathbf{y}=(y_{T-1},y_{T-2},\dots,y_{2},y_{1})$. We denote the model parameters as $\boldsymbol{\beta}=(a,c,\eta)$ and $b$. Then, the overall likelihood of the time series except for initial values is given by the following joint probability distribution:

[TABLE]

Since the time series is modeled as the AR model and the network structure of the time series is assumed to be asymmetric, the above likelihood function can be divided into the following product of two conditional distributions:

[TABLE]

Therefore, testing whether the variance $\Sigma_{x|y}$ decreases significantly, by including directed interactions, is equivalent to test whether the conditional likelihood $L_{x|y}$ increases significantly.

II.4 Parameter estimation

Log-likelihood functions of the observed time series are given by

[TABLE]

where we have omitted the constant terms not related to model parameters, and $\boldsymbol{\xi}=(\xi_{T-1},\xi_{T-2},\dots,\xi_{2})$. Though our interest is the likelihood function $L_{x|y}(\boldsymbol{\beta},b|\mathbf{x},\mathbf{y})$, this likelihood is over-parameterized and unidentifiable. Thus, we cannot determine the model parameters uniquely by using $L_{x|y}(\boldsymbol{\beta},b|\mathbf{x},\mathbf{y})$ alone. Therefore, we propose a parameter estimation procedure where we first estimate the parameter $b$ by using maximum likelihood estimation for $L_{y}(b|\mathbf{y})$ and then estimate the remaining parameters in $L_{x|y}(\boldsymbol{\beta},b|\mathbf{x},\mathbf{y})$ by substituting the previously estimated parameter $b=\hat{b}$. In other words, we first estimate the noise term $\xi_{t}$ by using maximum likelihood estimation for $L_{y}(b|\mathbf{y})$, then we substitute the maximum likelihood estimates $\xi_{t}=\hat{\xi}_{t}$ for the explanatory variables in $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$. The justification for this parameter estimation procedure is given as follows. Since we assume the asymmetric network structure, the model parameter in $L_{y}(b|\mathbf{y})$ can be estimated independently on $x_{t}$. That is, the model parameter $b$ can be estimated by $y_{t}$ alone, and we can obtain the residual time series $\hat{\xi}_{t}\equiv y_{t}-\hat{b}y_{t-1}$. Inserting the residual time series $\hat{\xi}_{t}$ into the likelihood function $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$ reduces the model for $x_{t}$ to a simple AR model with the explanatory variables $\hat{\xi}_{t}$ added, as shown in Eq. (27). Therefore, the maximum likelihood estimate of the inserted likelihood $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\hat{\boldsymbol{\xi}})$ always increases by including interaction term $c$.

Partial Granger causality in the literature suffers from a negative value of the test statistics. In Ref. Roelstraete and Rosseel (2012), it was pointed out that the reason for this may be due to small coding errors. We infer the reason for this from the discussion of the parameter estimation of the model. In the method of partial Granger causality, maximum likelihood estimation for joint distribution $p(x,z)$ was proposed as a method of parameter estimation Barrett et al. (2010), where $x$ and $z$ are the target variable and the conditioning variable, respectively. Even in our Granger causality measure, maximum likelihood estimation for the overall likelihood $L$ is a possible candidate for the parameter estimation. However, when we use maximum likelihood estimation for $L$, $L$ always increases by including the interaction term although $L_{x|y}$ does not necessarily increase. It is because, the interaction parameter included is used to increase $L$ itself, but not increase $L_{x|y}$. We have numerically confirmed that only occasionally maximum likelihood estimation for $L$ may lead to a negative value of the test statistics.

If the noise term $\xi_{t}$ were observed, the distribution of test statistics for inferring directed interactions can be obtained analytically based on the log-likelihood ratio statistic in the framework of the generalized linear model as follows. We denote the model parameters of our interest as $\boldsymbol{\beta}_{1}=(a,c,\eta)$, and the parameters of a full model as $\boldsymbol{\beta}_{\mathrm{max}}$. The full model is a generalized linear model with the same probability distribution and the same link function as the model of our interest, except that it has a maximum number of parameters that can be estimated from the data. The number of parameters for the model of our interest and the full model is $p_{1}$ and $T-1$, respectively. Since the likelihood function $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$ is expressed as a product of that of the Gaussian distribution, the deviance statistic $D_{1}$ obeys the chi-square distribution, $D_{1}\sim\chi^{2}(T-1-p_{1})$, provided that the model of our interest describes the data as much as the full model Dobson and Barnett (2008); Mccullagh and Nelder (2008):

[TABLE]

where

[TABLE]

and hats denote maximum likelihood estimates for $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$. On the other hand, if there is a significant difference in describing the data between model of our interest and the full model, the deviance obeys the non-central chi-square distribution and thus, takes a larger value than that expected by the central chi-square distribution.

When we want to test the existence of directed interactions, we set a null hypothesis $H_{0}$ and an alternative hypothesis $H_{1}$. Model parameters corresponding to these hypotheses are expressed as $\boldsymbol{\beta}_{0}=(a,\eta)$ and $\boldsymbol{\beta}_{1}=(a,c,\eta)$, respectively. These model are nested, that is, they have the same probability distribution and the same link function, but the parameters of the model for $H_{0}$ is a special case of the parameters of the model for $H_{1}$. The number of parameters for these models is $p_{0}$ and $p_{1}$, respectively. Then, we consider the difference of deviance between these models,

[TABLE]

where

[TABLE]

and hats denote maximum likelihood estimates for $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$. The difference of deviance obeys the central chi-square distribution, $\Delta D\sim\chi^{2}(p_{1}-p_{0})$, provided that both models describe the data as much as the full model. On the other hand, if there is a significant difference in describing the data between these models, the difference of deviance obeys the non-central chi-square distribution, and takes a larger value than that expected by the central chi-square distribution.

Here, we are assuming that the dispersion parameter of the model $\tau$ cannot be estimated from the data. Then, the difference of deviance is proportional to unknown parameter $\tau$, i.e., the difference of deviance is expressed as a scaled deviance. In practice, we can remove this unknown parameter in the test statistics by dividing deviances. First, we assume that the model for the alternative hypothesis describes the data as much as the full model. That is, the deviance $D_{1}$ obeys the central chi-square distribution. Next, we divide the difference of deviance $\Delta D$ by $D_{1}$ to make a new statistic that is independent of the unknown parameter $\tau$,

[TABLE]

The statistic $F$ obeys the $F$-distribution, $F\sim F(p_{1}-p_{0},T-1-p_{1})$, provided that both models describe the data as much as the full model. On the other hand, if there is a significant difference in describing the data between the null model and the alternative model, the statistic $F$ obeys the non-central $F$-distribution and takes a larger value than that expected by the central $F$-distribution. In summary, the existence of directed interactions can be tested as follows. If the test statistic $F$ is in the significance level of the $F$-distribution, we judge that there is no difference in describing the data between these models and accept a simpler model $H_{0}$. If the test statistic $F$ takes a value larger than the significance level, we reject $H_{0}$ and accept $H_{1}$.

The above result strictly holds true if $\xi_{t}$ were observed. In fact, $\xi_{t}$ is not observed and has to be estimated from the data. Here, we approximate $\xi_{t}$ by its maximum likelihood estimates for $L_{y}(b|\mathbf{y})$:

[TABLE]

Then, model parameters $\boldsymbol{\beta}$ are estimated from the approximated likelihood $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\hat{\boldsymbol{\xi}})$. In this paper, we assume that the test statistic obeys the $F$-distribution approximately even if we substitute $\hat{\xi}_{t}$.

The above result can be directly extended to a more general model where the model order of the interacting term and the lag of the noise correlation take any values. In this case, the model of the time series is expressed as follows:

[TABLE]

We express the model parameters as $\boldsymbol{\beta}=(\mathbf{a},\mathbf{c},\eta_{l_{\eta}})$, where $\mathbf{a}=(a_{1},a_{2},\cdots,a_{l_{a}})$, $\mathbf{c}=(c_{1},c_{2},\cdots,c_{l_{c}})$ and $l_{\eta}$ denotes the time point of lag for noise correlation and takes an integer value larger than zero. Then, the log-likelihood functions for this model are given by

[TABLE]

where

[TABLE]

It should be noted that since we are dealing with the AR model, the noise $\varepsilon_{t}$ and $\xi_{t}$ does not have autocorrelation. Therefore, $\varepsilon_{t}$ and $\xi_{t}$ cannot correlate with multiple time points. Thus, we just consider the correlation between $\varepsilon_{t}$ and $\xi_{t}$ with a single time point $\eta_{l_{\eta}}\xi_{t-l_{\eta}}$ rather than multiple time points $\sum_{i=1}^{l_{\eta}}\eta_{i}\xi_{t-i}$. The specific expression for the test statistic $F$ is given by

[TABLE]

where

[TABLE]

and hats denote maximum likelihood estimates for $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$. The statistic $F$ obeys the $F$-distribution, $F\sim F(p_{1}-p_{0},T-l_{\mathrm{max}}-p_{1})$, provided that both models describe the data as much as the full model.

Again, the model parameters $\mathbf{b}$ are calculated by maximum likelihood estimation for $L_{y}(\mathbf{b}|\mathbf{y})$ first, and the estimated values for $\xi_{t}$ is obtained as

[TABLE]

Then, the likelihood function $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$ is approximated by inserting $\xi_{t}=\hat{\xi}_{t}$:

[TABLE]

The model parameters $\boldsymbol{\beta}$ are estimated from the approximated likelihood $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\hat{\boldsymbol{\xi}})$, and we assume that the approximated test statistic approximately obeys the $F$-distribution, even if we substitute $\hat{\xi}_{t}$.

III Numerical experiments

In this section, we confirm the validity of our parameter estimation procedure described in the previous section and the distribution that the test statistic $F$ obeys numerically. Then, we propose a practical procedure for inferring directed interactions, which is based on a multiple statistical test for the case where the influence of latent common inputs may exist. In all numerical experiments in this paper, we assume that the model order of all the self interactions $l_{a}$ and $l_{b}$ is known for simplicity. In addition, we assume that the network structure is asymmetric, i.e., the interaction from $x_{t}$ to $y_{t}$ is absent.

We consider toy models which have been extensively applied in tests of Granger causality Baccalá and Sameshima (2001); Ding et al. (2006), except that we modify this model by adding latent common inputs to each time series. We simulated four toy models with and without latent common inputs and directed interactions, $(a)$, $(b)$, $(c)$ and $(d)$:

[TABLE]

where the noise has a following covariance matrix:

[TABLE]

The parameters $a_{1}$, $a_{2}$, $b_{1}$, $b_{2}$, $\sigma_{x}$ and $\sigma_{y}$ take the same value at each model, however, the other parameters take different values. The values of the model parameters are listed in TABLE 1. The schematic view of the simulated models and the run-sequence plots of time series are shown in FIG. 1 and FIG. 2, respectively.

III.1 Sampling distribution of the test statistics

In this subsection, we first confirm the influence of the correlated noise generated by latent common inputs numerically. For simplicity, we assume that the model order of directed interactions and the lag of correlated noise are known. For the case where these lags are unknown will be treated in the next subsection. We first confirm that the latent common inputs with delays produce a spurious directed interaction while the equal-time correlated noise does not, though only the latent common inputs without delays were considered in the literature. FIG. 3 shows the effect of delays of latent common inputs on the sampling distributions of the test statistics for normal Granger causality. The Figures show that when latent common inputs reach simultaneously, the test statistic is the same as that of the null hypothesis $F$-distribution, while in the case of latent common inputs with delay the test statistic deviates from the null distribution. The purpose of this paper is to construct a robust statistical measure for Granger causality against the latent common inputs with delays.

Next, we confirm the validity of our statistical measure and the parameter estimation procedure proposed in the previous section. FIG. 4 shows the test statistics for proposed Granger causality at model $(b)$ with various sample sizes. From the figure, the proposed test statistics approximately obey the $F$-distribution, and this property does not change even if the sample size varies. Thus, the proposed method is robust against the influence of latent common inputs. The reason for why the test statistic does not exactly obey the $F$-distribution is that we have used the approximated likelihood $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\hat{\boldsymbol{\xi}})$ rather than $L_{x|y}(\boldsymbol{\beta}|\mathbf{x},\mathbf{y},\boldsymbol{\xi})$.

The above result alone cannot deny the criticism that the proposed method is just insensitive to directed interactions. Therefore, we show that the proposed method has an ability to detect true directed interactions. FIG. 5 shows the test statistics of simulated model $(d)$ for various sample sizes. As the sample size increases, the test statistic deviates from the $F$-distribution, thus the method has an ability to detect directed interactions.

III.2 Practical procedure

The assumption made in the previous subsection that the model order of auto-regressive coefficients and the lag of correlated noise are known, is not realistic. In practice, these lags in true models are unknown in many cases. Furthermore, since the result of a statistical test depends crucially on the number of parameters to be tested, it is necessary to estimate these lags properly in order to infer the interactions accurately. Therefore, in this subsection, we propose a practical procedure for testing directed interactions, which is based on a stepwise variable increasing method through multiple tests.

We first estimate model parameter $y_{t}$ by maximum likelihood estimation for $L_{y}$, and obtain the residual time series $\hat{\xi}_{t}$. Then, we select the best lags of the noise correlation and the directed interaction, $l_{\eta}$ and $l_{c}$, by the Bayesian Information criterion (BIC) separately. BIC is given by

[TABLE]

where $l$ is the maximum number of lags to be searched for, $p$ is the number of model parameters and we have omitted the irrelevant constants which do not depend on the model parameters. Here we have assumed that when we calculate BIC, the dispersion parameter can be estimated from the data, $\tau=\hat{\tau}=\frac{1}{T-l}\sum_{t=l+1}^{T}(x_{t}-\hat{\mu}_{t})^{2}$. Next, using the selected lags which minimize BIC, the $p$-values for the noise correlation and the directed interaction are calculated separately. Then, if the $p$-value, which is the smaller one, takes a smaller value than the multiple test threshold, we reject the null hypothesis and adopt the alternative hypothesis. If either alternative hypothesis is adopted, we again select the lag and perform a statistical test for the other variable under the adopted alternative hypothesis. Since this procedure can be seen as multiple tests, we apply the Bonferroni correction to the threshold of the $p$-value in order for the familywise error rate to be less than the prescribed threshold.

We compare the performance in detecting directed interactions between the normal Granger causality procedure and the proposed procedure in FIG. 6. In this numerical experiment and below, we set the maximum number of lags to be searched, $l_{\eta}$ and $l_{c}$, be 6. Since the driving force of the AR model is the noise, the correlation between the variable $y_{t}$ and the noise $\xi_{t}$, in many cases, is large. Then, when we add both $y_{t}$ and $\xi_{t}$ on explanatory variables when exploring the best lags, the parameter estimation occasionally cannot perform due to multicollinearity. Therefore, if the multicollinearity occurs, we stop searching for lags. As in the case of the proposed procedure, we select the best number of lags for the directed interaction and perform a statistical test on it for the normal Granger causality procedure. We set the threshold of the $p$-value be 0.05 and the Bonferroni correction coefficient to $1/2$ in the case of the multiple test in order the familywise error rate to be 0.05.

From FIG. 6, it turns out that when the sample size is small, latent common inputs may produce a spurious directed interaction, even in the proposed method. It is because, since the noise is the driving force of the AR model, then both $p$-values for interaction $y_{t}$ and $\xi_{t}$ take small values, and a comparison of $p$-values between $y_{t}$ and $\xi_{t}$ fails. In addition, the proposed method has a lesser statistical power in detecting the true directed interaction than the normal Granger causality procedure. However, these shortcomings are resolved by increasing the sample size, and the directed interaction is detectable with high accuracy in any simulated models. On the other hand, the normal Granger causality procedure detects a spurious directed interaction with a high probability even in the small sample size, and inevitably detects a spurious directed interaction when the sample size is large.

In order to quantify the shortcoming of our proposed procedure in detecting the true directed interaction, we compare the statistical power in detecting directed interactions of the simulation model $(c)$ with the normal Granger causality procedure. FIG. 7 plots the statistical power against the sample size for both procedures. The proposed procedure has always poor statistical power than the normal Granger causality, i.e., the detection performance is conservative. However, the shortcoming is resolved by increasing the sample size, since the statistical power saturates to 1 as the sample size increases.

IV Conclusions

The original Granger causality test detects a spurious directed interaction when there are the confounding influence of latent common inputs and latent variables. In this paper, we propose a new Granger causality measure which is robust against the influence of latent common inputs. Our method is inspired by partial Granger causality in the literature and its variant, in which the influence of the latent common inputs between the source variable and the target variables was eliminated from the test statistics. We propose a statistical measure from which the correlated component of the noise between the source variable and the target variable was removed, which was caused by the influence of latent common inputs, while in partial Granger causality the correlated component of the noise between the target variable and the conditioning variable was removed. Then, we discuss the sampling distribution of our test statistics by using the likelihood function of time series and show numerically that the statistics approximately obey the $F$-distribution. According to what the authors in the literature say, the method of eliminating the correlated component of the noise between the target variable and the conditioning variable was called a partial Granger causality, since it was inspired by the partial correlation in statistics. However, in the light of the discussion in this paper, our method which eliminates the correlated component of the noise between the source variable and the target variable should be said to be a natural extension of Granger causality to the case where the noise has a correlation at a various time lag.

In practice, the model order of interactions and the lag of the correlated noise are unknown. Therefore, we propose a practical procedure in detecting directed interactions for the case where the influence of latent common inputs may exist, which is based on a stepwise variable increasing method by multiple tests. The performance of our procedure is verified by using numerically simulated time series. The normal Granger causality test inevitably detects a spurious directed interaction as the sample size increases, while the proposed method does not. The robustness of the proposed method against the influence of latent common inputs does not change as the sample size increases, thus our method can remove the influence of the latent common inputs properly. As a shortcoming of this robustness, the ability of proposed method in detecting the true directed interaction is weaker than the normal Granger causality procedure. However, this shortcoming can be resolved adequately by increasing the sample size, therefore the proposed method is practical enough.

Of the confounding environmental influence, this study only treats the latent common inputs while the influence of latent variables was not considered. It is because, the influence of latent variables appears as an auto-correlated noise, which cannot be dealt with in the framework of the AR model strictly. Granger causality in a model with an auto-correlated noise has to be defined in the framework of the MA model. The study in this direction has been worked on in the framework of the state space model Barnett and Seth (2015).

Our method assumes that the network structure of the interaction between each variable is asymmetric. Therefore, the extension to the case where the networks structure may be symmetric will be in future works. In addition, it is important to show a theoretical justification of the parameter estimation proposed in this paper. Furthermore, it is interesting to extend the proposed method to point process models, such as spike activity model of neuron.