{{ message }}

Modeling: Granger causality analysis

Fa-Hsuan Lin edited this page Nov 2, 2019 · 8 revisions
Clone this wiki locally

Granger causality is the method of inferring directional interaction between time series, which in the context of neuroimaging can be derived from statistical parametric mapping of fMRI/MEG/EEG measurements. Most Granger causality values are calculated using auto-regressive (AR) modeling of the time series. Here are some examples in Matlab.

The multivariate auto-regressive modeling toolbox (ARFIT) is very useful and has been extensively applied in many studies. Our codes are built on top of this toolbox.

Here are plots of the dynamic Granger causality estimates for propagated spikes.

Codes

The function etc_granger.m calculates the Granger causality values and the associated statistics

Examples

Simulating bi-variate AR model time series

Here we generate bi-variate time series and calculate the Granger causality. This is the example of Granger causality values from bi-variate time series generated from a 1-st order AR model in Roebroeck et al. (NeuroImage 2005) without convolving the time series with hemodynamic response function. Mathematically, it is

``````[x(t+1), y(t+1)]^T=[-0.9, 0; I, -0.9][x(t), y(t)]^T+[ex(t), ey(t)], where [ex(t) ey(t)]^T are bi-variate noise with a covariance matrix S=[1 0; 0 1];
``````

I is the coupling strength and it is set I=0.5 for strong influence. The model clearly suggests a causal modulation from x -> y.

code: granger_ar_sim.m

Simulating the bi-variate AR model time series

Here we generate dynamic bi-variate time series and calculate the dynamic Granger causality. This is the example of Granger causality values from bi-variate time series generated from a 2nd order AR model. Mathematically, it is

``````[x(t+2), y(t+2)]^T=[-0.9, 0; I, -0.9][x(t+1), y(t+1)]^T+[-0.3, 0; I*0.5, -0.3][x(t), y(t)]^T+[ex(t), ey(t)], where [ex(t) ey(t)]^T are bi-variate noise with a covariance matrix S=[1 0; 0 1];
``````

I is the coupling strength and it is set I=0.5 for strong influence. Such bi-variate time series is only active between time sample 300 and 600 among total 1000 samples. Other time instants are two white noises. Specifically, the final time series [X(t), Y(t)]^T is [X(t), Y(t)]^T = env(t).[x(t), n(t)]^T + (1-env(t)).[xn(t), yn(t)]^T, where env(t) is a time course with value close to 1 between sample 300 and 600 (sigmoidal ramping) and 0 else where. [xn(t) yn(t)]^T are white Gaussian noise with unit variance.

The goal here is to estimate the appropriate window length and AR model order using the SURE metric. Subsequently, dynamic Granger causality is estimated using the chosen AR model order and window length.

code: dgranger_ar_sim.m