Granger causality analysis with nonuniform sampling and its application to pulse-coupled nonlinear dynamics

The Granger causality (GC) analysis is an effective approach to infer causal relations for time series. However, for data obtained by uniform sampling (i.e., with an equal sampling time interval), it is known that GC can yield unreliable causal inference due to aliasing if the sampling rate is not sufficiently high. To solve this unreliability issue, we consider the nonuniform sampling scheme as it can mitigate against aliasing. By developing an unbiased estimation of power spectral density of nonuniformly sampled time series, we establish a framework of spectrum-based nonparametric GC analysis. Applying this framework to a general class of pulse-coupled nonlinear networks and utilizing some particular spectral structure possessed by these nonlinear network data, we demonstrate that, for such nonlinear networks with nonuniformly sampled data, reliable GC inference can be achieved at a low nonuniform mean sampling rate at which the traditional uniform sampling GC may lead to spurious causal inference.

Figure 1

The GC sampling structure (GC value as a function of sampling interval length $\tau$). The GC analysis is applied to the time series obtained from a two-neuron I&F network with a unidirectional connection from neuron $y$ to neuron $x$. Plotted are $F_{x\to y}$ (red) and $F_{y\to x}$(cyan) estimated from voltage time series of the system (1) through linear regression. Here, $\mu =1\mathrm{kHz},\lambda =0.03$, with coupling strength $s=0.02$ ($s_{ij}=s$ for $s_{ij}\ne 0$). (a) The oscillatory behavior of GC vs $\tau$. Shaded regions represent a $95\%$ confidence interval. (b) The vanishing of GC as $\tau$ approaches 0. Note that the estimation biases of GC are removed in (b) (see Refs. [20, 21] for details).

Figure 2

The dependence of GC on the cutoff frequency $f_{\mathrm{cut}}$ for nonuniform sampling. (a) GC vs $f_{\mathrm{cut}}$ and GC vs $1/2f_{\mathrm{cut}}$ (inset), (b) normalized GC vs $f_{\mathrm{cut}}$ and normalized GC vs $1/2f_{\mathrm{cut}}$ (inset). Plotted are (a) $F_{y\to x}$ (cyan), $F_{x\to y}$ (red), and (b) $2f_{\mathrm{cut}}{F}_{y\to x}$ (cyan), $2f_{\mathrm{cut}}{F}_{x\to y}$ (red) as a function of $f_{\mathrm{cut}}$. (c) Coherence as a function of $f$. For Figs. 2, 3, 4, 5, 6, the time series are generated from the same two-neuron I&F network with a unidirectional connection from neuron $y$ to neuron $x$ with parameters $\mu =1\mathrm{kHz},\lambda =0.012$, and coupling strength $s=0.02$ ($s_{ij}=s$ for $s_{ij}\ne 0$). Note that here we use PSD functions that are accurately evaluated from a sufficiently long time series (${10}^8$ ms) and sufficiently small sampling interval length $\left(0.0625\mathrm{ms}\right)$ for GC analysis. These PSD functions are regarded as the empirical PSD. The black arrows in (a)–(c) point to $f_{\mathrm{cut}}=1\mathrm{kHz}$.

Figure 3

Comparison between estimated PSDs: empirical PSDs (red), nonuniform sampling PSDs (black), nonuniform sampling PSDs after covariance truncation (cyan), nonuniform sampling PSDs after covariance truncation, and high-frequency power-law tail-fitting (dotted blue). (a) The auto-PSD for neuron y: $P_{yy}$. (b) The real part of the cross PSD: $\mathrm{Re}\left(P_{yx}\right)$. (c) The imaginary part of the cross PSD: $\mathrm{Im}\left(P_{yx}\right)$. Insets are the corresponding log-log plots of their absolute values. The PSDs are estimated from $X_t$ and $Y_t$ generated by a two-neuron I&F network with a unidirectional connection from neuron $y$ to neuron $x$.

Figure 4

Exponential decay of covariances. $\left|C_{xx}\right|$ (red), $\left|C_{xy}\right|$ (black), $\left|C_{yy}\right|$ (cyan) are estimated from bivariate time series $X_t,Y_t$ generated by a two-neuron I&F model with a unidirectional connection from neuron $y$ to neuron $x$.

Figure 5

Covariances and their derivatives. (a) Auto-covariance $C_{xx}\left(\tau \right)$, (b) $\frac{d{C}_{xx}\left(\tau \right)}{d\tau }$, (c) cross-covariance $C_{xy}\left(\tau \right)$, (d) $\frac{d^2{C}_{xy}\left(\tau \right)}{d{\tau }^2}$, computed from the voltage time series $X\left(t\right)$ and $Y\left(t\right)$ of neuron $x$ and neuron $y$, respectively, which are generated by a two-neuron I&F network with a unidirectional connection from neuron $y$ to neuron $x$. From the full simulation of the I&F dynamics of neuron $x$, we obtain the firing rate $R_x\approx 0.065m{s}^{-1}$. From numerical computation of the auto-covariance $C_{xx}\left(\tau \right)$, we have ${\left.\frac{d{C}_{xx}\left(\tau \right)}{d\tau }\right|}_{0^{+}}=-{\left.\frac{d{C}_{xx}\left(\tau \right)}{d\tau }\right|}_{0^{-}}=-0.032{\mathrm{ms}}^{-1}$, which is equal to $-\frac{1}{2}{\left(x^{\mathrm{th}}-x^{\mathrm{r}}\right)}^2{R}_x$ within our numerical accuracy.

Figure 6

Comparison between the normalized GC sampling structures. $F_{x\to y}$ (dash), $F_{y\to x}$ (solid) are obtained by the conventional linear regression (cyan), by the uniform sampling nonparametric GC analysis (blue), by the nonuniform sampling GC analysis with truncation and fitting procedures (red). The time series are generated by a two-neuron I&F network with a unidirectional connection from neuron $y$ to neuron $x$. Note that, for the evaluation of each normalized GC point in the figure, we use exactly the same number of data points of the time series for every method.

Figure 7

Conditional nonuniform GC analysis for a 10-neuron I&F network. The parameters are $\mu =1\mathrm{kHz},\lambda =0.01$, with coupling strength $s=0.01$($s_{ij}=s$ for $s_{ij}\ne 0$). (a) The adjacency matrix $A_{ij}$ of the network. $A_{ij}=0$ (white) when $s_{ij}=0$, and $A_{ij}=1$ (black) when $s_{ij}\ne 0$. (b) Ranked normalized GC values obtained through conditional nonuniform GC framework with truncation and fitting procedures. The nonuniform mean sampling interval length is $\tau =10\mathrm{ms}$. Blue circles indicate the GC values for uncoupled directions whereas red circles indicate the GC values for coupled directions. The solid black line indicates a natural GC threshold $\left(8.5\times {10}^{-5}\right)$, below which the corresponding direction is inferred as uncoupled. By applying such a GC threshold, we can well reconstruct the network topology of the I&F system.