Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improvements to Granger Causality evaluation #428

Closed
wants to merge 9 commits into from
+258 −304
Closed

Improvements to Granger Causality evaluation #428

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
52f5cdf
Improvements to Granger Causality evaluation
DavideNuzzi Jul 25, 2021
b5d3f1e
Delete YuleWalker_Mask.m
DavideNuzzi Aug 18, 2021
131c1e2
Delete YuleWalker_Inverse.m
DavideNuzzi Aug 18, 2021
17b3c57
Delete bst_granger_conditional.m
DavideNuzzi Aug 18, 2021
727b048
Delete bst_granger.m
DavideNuzzi Aug 18, 2021
c98de7d
Delete bst_granger_spectral.m
DavideNuzzi Aug 18, 2021
7373c86
Delete bst_granger_spectral_conditional.m
DavideNuzzi Aug 18, 2021
e0014f3
Add files via upload
DavideNuzzi Aug 18, 2021
1e7ca9e
Add files via upload
DavideNuzzi Aug 18, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
246 changes: 63 additions & 183 deletions toolbox/connectivity/bst_granger.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
function [connectivity, pValue, connectivityV, pValueV, X, Y] = bst_granger(X, Y, order, inputs)
% BST_GRANGER Granger causality in mean and variance between any two
% signals, using two Wald statistics
% in mean: regular log-GC from Geweke1982
% in variance: information statistic from Hafner2007
function [connectivity, pValue] = bst_granger(X, Y, order, inputs)
% BST_GRANGER Granger causality between any two signals,
% using two Wald statistics
%
% Inputs:
% sinks - first set of signals, one signal per row
Expand All @@ -20,58 +18,32 @@
% |-flagFPE - if true, optimize order for AR model
% | if false (default), force same order in all AR models
% | [E: default false]
% |-lag - maximum lag in ARCH model for causality in variance
% | [S: nonnegative integer]
% |-flagELM - if true, optimize order for ARCH model
% | if false (default), force same order in all ARCH models
% | [L: default false]
% |-rho - ADMM parameter from augmented Lagrangian
% | --> lower means faster but at cost of stability
% | --> higher means convergence but at cost of speed
% | --> 50 is a good starting point
% | [R: nonnegative number, default = 50]
%
% Outputs:
% connectivity - A x B matrix of causalities in mean from source to sink
% connectivity - A x B matrix of causalities from source to sink
% [C: MX x MY matrix]
% pValue - parametric p-value for corresponding Granger causality in
% mean estimate
% pValue - parametric p-value for corresponding Granger causality estimate
% [P: MX x MY matrix]
% connectivityV - A x B matrix of causalities in variance from source to sink
% [CV: MX x MY matrix]
% pValueV - parametric p-value for corresponding Granger causality in
% variance estimate
% [PV: MX x MY matrix]
%
% See also BST_MVAR, BST_VGARCH.
%
% For each signal pair (a,b) we calculate the Granger causality in mean GC(a,b):
% Var(x_a[t] | x_a[t-1, ..., t-k])
% ----------------------------------------------------
% Var(x_a[t] | x_a[t-1, ..., t-k], y_b[t-1, ..., t-k])
% If Y is empty or Y = X, we set element GC(a,a) to be zero.
%
% If inputs.lag is set, then for each signal pair (a,b) we calculate the Wald
% statistic EC(a,b) testing whether C_{a,b} = 0 where
% x[n] = A[1] x[n-1] + ... + A[P] x[n-P] + e[n]
% e[n] ~ normal with mean 0 and covariance H[n]
% H[n] = W*W' + sum_{r=1}^{inputs.lag} C_r' * e[n-r] * e[n-r]' * C_r
% If Y is empty or Y = X, we set element EC(a,a) to be zero.
% See also BST_MVAR
%
% For each signal pair (a,b) we calculate the Granger causality:
%
% det(restricted_residuals_cov_matrix)
% gc = ----------------------------------------------------
% det(full_cov_matrix)
%
% see Cohen, Dror, et al. "A general spectral decomposition of causal influences
% applied to integrated information." and Barret, Barnett and Seth "Multivariate
% Granger causality and generalized variance".
%
% Call:
% [inMean, inVariance] = bst_granger(X, Y, 5, inputs);
% [inMean, inVariance] = bst_granger(X, [], 5, inputs); % every pair in X
% [inMean, inVariance] = bst_granger(X, [], 20, inputs); % more delays in AR
% connectivity = bst_granger(X, Y, 5, inputs);
% connectivity = bst_granger(X, [], 5, inputs); % every pair in X
% connectivity = bst_granger(X, [], 20, inputs); % more delays in AR
% inputs.nTrials = 10; % use trial-averaged covariances in AR estimation
% inputs.standardize = true; % zero mean and unit variance
% inputs.flagFPE = true; % allow different orders for each pair of signals
% inputs.lag = 3; % estimate causality in variance using lag-3 ARCH model
% inputs.flagELM = true; % find the optimal lag rather than always given lag
% inputs.rho = 50; % parameter to tune estimation of ARCH model
% % higher -> more stability, much slower
% % lower -> might not converge, much faster
% % go in multiples of 5 up or down as desired
% % 50 is a good start for rho

% @=============================================================================
% This function is part of the Brainstorm software:
Expand All @@ -92,6 +64,7 @@
% =============================================================================@
%
% Authors: Sergul Aydore & Syed Ashrafulla, 2012
% Modified by: Davide Nuzzi, 2021

% default: 1 trial
if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials)
Expand All @@ -103,16 +76,6 @@
inputs.flagFPE = false;
end

% default: do not optimize order in MVAR modeling
if ~isfield(inputs, 'flagELM') || isempty(inputs.flagELM)
inputs.flagELM = false;
end

% default: ADMM works well with rho = 50
if ~isfield(inputs, 'rho') || isempty(inputs.rho)
inputs.rho = 50;
end

% dimensions of the signals
nX = size(X, 1);
if ndims(X) == 3
Expand Down Expand Up @@ -153,75 +116,48 @@

%% Iterate over all pairs of sinks & sources

% for causality in mean we need the restricted variance
restOrder = zeros(nX, 1); restCovFull = zeros(nX, order+1);
for iX = 1:nX
[syed, syed, restOrder(iX), syed, syed, restCovFull(iX, :)] = bst_mvar(X(iX, :), order, inputs.nTrials, inputs.flagFPE); %#ok<ASGLU>
end

if isempty(Y) % auto-causality

% setup
connectivity = zeros(nX);
connectivityV = zeros(nX);

% only iterate over one triangle
for iX = 1:nX

% iterate over all the pairs after iX
for iY = (iX+1):nX

% bivariate autoregressive model with sink_a and sink_b
[syed, syed, unOrder, syed, syed, unCovFull, residual] = bst_mvar([X(iX, :); X(iY, :)], order, inputs.nTrials, inputs.flagFPE); %#ok<ASGLU>

% causality in mean: Geweke-Granger, i.e. restricted variance / unrestricted variance - 1
if inputs.flagFPE % get the minimum order of the two models estimated

% source = iY, sink = iX
minOrder = min([restOrder(iX) unOrder]);
connectivity(iX, iY) = restCovFull(iX, minOrder+1) / unCovFull(1, 1, minOrder+1) - 1;

% source = iX, sink = iY
minOrder = min([restOrder(iY) unOrder]);
connectivity(iY, iX) = restCovFull(iY, minOrder+1) / unCovFull(2, 2, minOrder+1) - 1;

else % by default, bst_mvar sends the result of the single model of given order into the "Full" variables
% two-variate model
[transfers, noiseCovariance, order] = bst_mvar([X(iX, :); X(iY, :)], order, inputs.nTrials, inputs.flagFPE);

% data correlations using Yule-Walker (up to high order 50)
R = yule_walker_inverse(transfers, noiseCovariance, 50);

% restricted model iY -> iX

connectivity(iX, iY) = restCovFull(iX) / unCovFull(1, 1) - 1; % source = iY, sink = iX
connectivity(iY, iX) = restCovFull(iY) / unCovFull(2, 2) - 1; % source = iX, sink = iY
% mask for the coefficients of the restricted model
mask = ones(2);
mask(1,2) = 0;

end
% restricted bivariate model (using masked row-by-row solution of
% YW equations)
[transfers_restricted,noiseCovariance_restricted] = yule_walker_mask(R, mask);

% causality in variance
if isfield(inputs, 'lag') && inputs.lag > 0 && any(abs(residual(1, :) - residual(2, :)) > eps)

% preprocess the residual first
residual = bst_bsxfun(@rdivide, residual, std(residual, [], 2));

% bivariate ARCH estimation
[W, C, D, R, S, information, rhoBest] = ... % bivariate ARCH modeling
bst_vgarch('vec', residual, inputs.nTrials, inputs.lag, 0, inputs.flagELM, 3, 'on', 999, [], inputs.rho, false, [], [], [], []); %#ok<ASGLU>

% use SQP if ADMM failed (found by an exploding augmentation parameter rho)
if rhoBest > 10 * inputs.rho || any(isinf(information(:)))
try % SQP may fail too because the data is nonstationary; rho will be -1 if it succeeds, indicating ADMM did not work
[W, C, D, R, S, information] = ... % bivariate ARCH modeling
bst_vgarch('vec', residual, inputs.nTrials, inputs.lag, 0, inputs.flagELM, 2, 'off', 999, [], [], false, [], [], [], []); %#ok<ASGLU>
catch ME %#ok<NASGU> % this data kills ADMM & SQP so we assume no causality in variance
W = bst_correlation(residual, [], struct('normalize', false, 'nTrials', 1, 'maxDelay', 0, 'nDelay', 1, 'flagStatistics', false)); % W is covariance
C = zeros(2*(2+1)/2, 2*(2+1)/2*inputs.lag); information = eye(2*2 + 2*(2+1)/2*2*(2+1)/2*inputs.lag); % and all the vARCH parameters are zero
end
end

% Wald statistic
theta = [W(:); C(:); D(:)]; % stack into parameter vector
covariance = inv(information); % estimated covariance matrix
J = sort([7 + 9*(0:inputs.lag-1) 10 + 9*(0:inputs.lag-1)]); % indices corresponding to elements C_{r,12} and C_{r,13} for all lags r in vARCH model
connectivityV(iX, iY) = theta(J)' / covariance( J, J) * theta(J); % source -> sink: use C_{12}=7, C_{13}=10 & add 9 for the other lags
connectivityV(iY, iX) = theta(J-1)' / covariance(J-1,J-1) * theta(J-1); % sink -> source: use C_{32}=9, C_{31}=6 & add 9 for the other lags
% connectivity
connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance));

% restricted model iX -> iY

end
% mask for the coefficients of the restricted model
mask = ones(2);
mask(2,1) = 0;

% restricted bivariate model (using masked row-by-row solution of
% YW equations)
[transfers_restricted,noiseCovariance_restricted] = yule_walker_mask(R, mask);

% connectivity
connectivity(iY, iX) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance));
end

% diagonal will equal the maximum of all inflows and outflows for iX
Expand All @@ -233,53 +169,28 @@

% setup
connectivity = zeros(nX, nY);
connectivityV = zeros(nX, nY);
duplicates = zeros(0, 2);

for iX = 1:nX
for iY = 1:nY

if any(abs(X(iX, :) - Y(iY, :)) > eps) % avoid duplicates

% bivariate autoregressive model with sink_a and source_b
[syed, syed, unOrder, syed, syed, unCovFull, residual] = bst_mvar([X(iX, :); Y(iY, :)], order, inputs.nTrials, inputs.flagFPE); %#ok<ASGLU>

% causality in mean: Geweke-Granger, i.e. restricted variance / unrestricted variance - 1
if inputs.flagFPE % get the minimum order of the two models estimated
minOrder = min([restOrder(iX) unOrder]);
connectivity(iX, iY) = restCovFull(iX, minOrder+1) / unCovFull(1, 1, minOrder+1) - 1;
else % by default, bst_mvar sends the result of the single model of given order into the "Full" variables
connectivity(iX, iY) = restCovFull(iX) / unCovFull(1, 1) - 1;
end

% causality in variance
if isfield(inputs, 'lag') && inputs.lag > 0 && any(abs(residual(1, :) - residual(2, :)) > eps)

% preprocess the residual first
residual = bst_bsxfun(@rdivide, residual, std(residual, [], 2));

% bivariate ARCH estimation
[W, C, D, R, S, information, rhoBest] = ... % bivariate ARCH modeling
bst_vgarch('vec', residual, inputs.nTrials, inputs.lag, 0, inputs.flagELM, 3, 'on', 999, [], inputs.rho, false, [], [], [], []); %#ok<ASGLU>

% use SQP if ADMM failed (found by an exploding augmentation parameter rho)
if rhoBest > 10 * inputs.rho || any(isinf(information(:)))
try % SQP may fail too because the data is nonstationary; rho will be -1 if it succeeds, indicating ADMM did not work
[W, C, D, R, S, information] = ... % bivariate ARCH modeling
bst_vgarch('vec', residual, inputs.nTrials, inputs.lag, 0, inputs.flagELM, 2, 'off', 999, [], [], false, [], [], [], []); %#ok<ASGLU>
catch ME %#ok<NASGU> % this data kills ADMM & SQP so we assume no causality in variance
W = bst_correlation(residual, [], struct('normalize', false, 'nTrials', 1, 'maxDelay', 0, 'nDelay', 1, 'flagStatistics', false)); % W is covariance
C = zeros(2*(2+1)/2, 2*(2+1)/2*inputs.lag); information = eye(2*2 + 2*(2+1)/2*2*(2+1)/2*inputs.lag); % and all the vARCH parameters are zero
end
end

% Wald statistic
theta = [W(:); C(:); D(:)]; % stack into parameter vector
covariance = inv(information); % estimated covariance matrix
J = sort([7 + 9*(0:R-1) 10 + 9*(0:R-1)]); % indices corresponding to elements C_{r,12} and C_{r,13} for all lags r in vARCH model
connectivityV(iX, iY) = theta(J)' / covariance( J, J) * theta(J); % source -> sink: use C_{12}=7, C_{13}=10 & add 9 for the other lags

end
% two-variate model
[transfers, noiseCovariance, order] = bst_mvar([X(iX, :); Y(iY, :)], order, inputs.nTrials, inputs.flagFPE);

% data correlations using Yule-Walker (up to high order 50)
R = yule_walker_inverse(transfers, noiseCovariance, 50);

% mask for the coefficients of the restricted model
mask = ones(2);
mask(1,2) = 0;

% restricted bivariate model (using masked row-by-row solution of YW equations)
[transfers_restricted,noiseCovariance_restricted] = yule_walker_mask(R, mask);

% connectivity
connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance));

else % save duplicates to modify later

Expand All @@ -300,39 +211,8 @@

%% Statistics: parametric p-values for causality in mean (based on regression coefficients) and variance (based on Wald statistics)

% causality in mean: F statistic of connectivity when multiplied by number of regressors
% F statistic of connectivity when multiplied by number of regressors
pValue = 1 - betainc(connectivity ./ (1 + connectivity), order / 2, (nSamples - order * inputs.nTrials - 2 * order - 1) / 2, 'lower');
% here we assume we have many more samples than the order of the MVAR model so that in all cases we use the second condition below to compute the p-value
% note: if connectivity = 0 (auto-causality or two of the same signals) then this formula evalutes pValue = 1 which is desired (no significant causality)

% causality in mean: F statistic of connectivity when multiplied by number of regressors
% tic
% iFlip = nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1 > connectivity * order;
% pValue = zeros(size(connectivity));
% pValue(~iFlip) = 1 - betainc(1 ./ (1 + connectivity(~iFlip)), (nSamples - order(~iFlip) * inputs.nTrials - 2 * order(~iFlip) - 1) / 2, order(~iFlip) / 2, 'upper');
% pValue(iFlip) = 1 - betainc(connectivity(iFlip) ./ (1 + connectivity(iFlip)), order(iFlip) / 2, (nSamples - order(iFlip) * inputs.nTrials - 2 * order(iFlip) - 1) / 2, 'lower');
% toc
% % fcdf(connectivity .* (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) ./ order, order, nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1)
% % which for nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1 <= connectivity * order is
% % = betainc((nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) ./ (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1 + connectivity * (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / order * order), (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / 2, order/2, 'upper')
% % = betainc((nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) ./ ((nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) .* (1 + connectivity)), (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / 2, order/2, 'upper')
% % = betainc(1 ./ (1 + connectivity), (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / 2, order/2, 'upper')
% % and for nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1 >= connectivity * order is
% % = betainc(connectivity .* (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) ./ order .* order ./ (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1 + connectivity .* (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) ./ order .* order), order/2, (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / 2, 'lower')
% % = betainc(connectivity .* (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) ./ ((nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) .* (1 + connectivity)), order/2, (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / 2, 'lower')
% % = betainc(connectivity ./ (1 + connectivity), order/2, (nSamples - order(~iFlip) * inputs.nTrials - 2 * order - 1) / 2, 'lower')

% causality in mean: chi-square statistic of connectivity when multiplied by the number of samples
% tic
% pValue = 1 - gammainc(connectivity * (nSamples - order * inputs.nTrials) / 2, order / 2);
% toc
% % chi2cdf(connectivity * (nSamples - order * inputs.nTrials), order)
% % = gamcdf(connectivity * (nSamples - order * inputs.nTrials), order/2, 2)
% % = gammainc(connectivity * (nSamples - order * inputs.nTrials) / 2, order / 2)

% causality in variance: chi-square statistic of connectivity when multiplied by number of samples (minus lag minus 1) to get chi-square statistic
pValueV = 1 - gammainc(connectivityV .* (nSamples - order * inputs.nTrials - inputs.lag - 1) / 2, inputs.lag);
% chi2cdf(connectivityV .* (nSamples - order * inputs.nTrials - inputs.lag - 1), 2 * inputs.lag)
% = gamcdf(connectivityV .* (nSamples - order * inputs.nTrials - inputs.lag - 1), (2 * inputs.lag)/2 = inputs.lag, 2)
% = gammainc(connectivityV .* (nSamples - order * inputs.nTrials - inputs.lag - 1) / 2, inputs.lag)
% note: if connectivityV = 0 (auto-causality or two of the same signals) then this formula evalutes pValueV = 1 which is desired (no significant causality)

Loading