title | tags | ||||
---|---|---|---|---|---|
Conditional Granger causality in the frequency domain |
|
Conditional Granger causality is a derivative of spectral Granger causality that is computed over a triplet of channels (or blocks of channels). It provides the advantage that for this triplet, it allows to differentiate between a delayed parallel drive from sources A to be B and C and a sequential drive from A to B to C.
This example illustrates the simulation and base analysis of the paper Chen, Y., Bressler, S. L., & Ding, M. (2006). Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. Journal of neuroscience methods, 150(2), 228-237.
See also: [the connectivity tutorial]/tutorial/connectivity/).
First, define parameters under which samples should be simulated.
simcfg = [];
simcfg.ntrials = 500;
simcfg.triallength = 1;
simcfg.fsample = 200;
simcfg.nsignal = 3;
simcfg.method = 'ar';
We want to simulate a system with three signals. Their noise is modeled as white noise processes with zero mean and standard deviations
σ12=1
σ22=0.2
σ32=0.3.
They will have the covariances ζ, η and ε. We also require a paramters μ=0.5.
parmeters of the model itself
mu = 0.5;
absnoise = [ 1.0 0.2 0.3 ];
First, we generate the sample for the case of sequential driving. We want to incorporate the system
x(t) = ζ(t)
y(t) = x(t-1) + η(t)
z(t) = μ⋅z(t-1) + y(t-1) + ε(t),
which we can do like this:
params(i,j,k): j -> i at t=k
simcfg.params(:,:,1) = [ 0 0 0;
1.0 0 0;
0 1.0 mu];
Note that the matrix representation for the covariance reads from columns to row, other than the MVAR-model is read intuitively. But we still need to hand the parameters of the noise to the model:
paper defines stds, not cov:
simcfg.noisecov = diag(absnoise.^2);
data2 = ft_connectivitysimulation(simcfg);
Now create sample data for the case of differentially delayed driving,
x(t) = ζ(t)
y(t) = x(t-1) + η(t)
z(t) = μ⋅z(t-1) + x(t-2) + ε,
which we can write as
simcfg.params(:,:,1) = [ 0 0 0;
1.0 0 0;
0 0 mu];
simcfg.params(:,:,2) = [ 0 0 0;
0 0 0;
1.0 0 0];
We build the actual MVAR-representation...
data1 = ft_connectivitysimulation(simcfg);
... and have a first look at the data:
figure
plot(data1.time{1}, data1.trial{1})
legend(data1.label)
xlabel('time (s)')
Don't be confused that we started with data2 and conclude with data1. This is just to maintain the order the systems have in the paper.
We generate spectral representations from the MVAR representations we defined with data1 and data2. After all, we want to compute spectral Granger causality. Fast Fourier is a good starting point.
freq = [];
freq.freqcfg = [];
freq.freqcfg.method = 'mtmfft';
freq.freqcfg.output = 'fourier';
freq.freqcfg.tapsmofrq = 2;
freqdata1 = ft_freqanalysis(freq.freqcfg, data1);
freqdata2 = ft_freqanalysis(freq.freqcfg, data2);
Let first compute regular bivariate Granger causality, as this makes the difference clear to what we want.
grangercfg = [];
grangercfg.method = 'granger';
grangercfg.granger.conditional = 'no';
grangercfg.granger.sfmethod = 'bivariate';
gdata = [];
gdata.g1_bivar_reg = ft_connectivityanalysis(grangercfg, freqdata1);
gdata.g2_bivar_reg = ft_connectivityanalysis(grangercfg, freqdata2);
However, we clearly want a multivariate approach. Also, we need to define channel combinations, as we now require triplets of inputs.
grangercfg.granger.conditional = 'yes';
grangercfg.channelcmb = {'signal001', 'signal002', 'signal003'};
grangercfg.granger.sfmethod = 'multivariate';
grangercfg.granger.conditional = 'yes';
block-wise causality
grangercfg.granger.block(1).name = freqdata1.label{1};
grangercfg.granger.block(1).label = freqdata1.label(1);
grangercfg.granger.block(2).name = freqdata1.label{2};
grangercfg.granger.block(2).label = freqdata1.label(2);
grangercfg.granger.block(3).name = freqdata1.label{3};
grangercfg.granger.block(3).label = freqdata1.label(3);
gdata.g1_multi_reg_conditional = ft_connectivityanalysis(grangercfg, freqdata1);
gdata.g2_multi_reg_conditional = ft_connectivityanalysis(grangercfg, freqdata2);
The label combinations are 6x2 cell arrays, containing all 2-permutations tuplets from the channels. How to interpret this? Is the combination a, b representing Fa→b|c? Let's check this. In scenario 2, we should clearly see a higher causality from 1→3 | 2 than in scenario 1 of the differentially delayed drive. This corresponds to row 4 in the gdata.g1_multi_reg_conditional.labelcmb. So, let's compare the labelcmb 1, 3 in both scenarios:
scenario1_mean = mean(gdata.g1_multi_reg_conditional.grangerspctrm(4, :));
scenario2_mean = mean(gdata.g2_multi_reg_conditional.grangerspctrm(4, :));