From 52f5cdfe20426b730b785470681f4870a8beb936 Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Sun, 25 Jul 2021 18:52:28 +0200 Subject: [PATCH 1/9] Improvements to Granger Causality evaluation Modified bst_granger.m and bst_granger_spectral.m (their results were not correct) Implemented the conditional GC both in time and frequency domain. --- toolbox/connectivity/YuleWalker_Inverse.m | 34 ++ toolbox/connectivity/YuleWalker_Mask.m | 38 ++ toolbox/connectivity/bst_granger.m | 246 ++++-------- .../connectivity/bst_granger_conditional.m | 236 ++++++++++++ toolbox/connectivity/bst_granger_spectral.m | 221 +++++------ .../bst_granger_spectral_conditional.m | 350 ++++++++++++++++++ 6 files changed, 821 insertions(+), 304 deletions(-) create mode 100644 toolbox/connectivity/YuleWalker_Inverse.m create mode 100644 toolbox/connectivity/YuleWalker_Mask.m create mode 100644 toolbox/connectivity/bst_granger_conditional.m create mode 100644 toolbox/connectivity/bst_granger_spectral_conditional.m diff --git a/toolbox/connectivity/YuleWalker_Inverse.m b/toolbox/connectivity/YuleWalker_Inverse.m new file mode 100644 index 000000000..48def697f --- /dev/null +++ b/toolbox/connectivity/YuleWalker_Inverse.m @@ -0,0 +1,34 @@ +function R = YuleWalker_Inverse(A,Sigma, maxOrder) +%YULEWALKER Summary of this function goes here +% A partire dalla matrice A e dalle correlazioni dei residui ottengo le +% correlazioni a vari ritardi. Le Gamma sono matrici N x N e in totale +% saranno (P+1) +% L'approccio dovrebbe essere di trasformare il problema in un'equazione +% di Lyapunov + + [N,~] = size(Sigma); + p = size(A,2)/N; + + % Copiato da CElinVAR + Im = eye(N*p); + F = [A; Im(1:end-N,:)]; + Delta=zeros(N*p); + Delta(1:N,1:N) = Sigma; + BigSigma = dlyap(F,Delta); + + R=NaN*ones(N,N,maxOrder+1); + for i=1:p + R(:,:,i)=BigSigma(1:N,N*(i-1)+1:N*i); + end + + for k=p+1:maxOrder+1 + Rk=R(:,:,k-1:-1:k-p); + Rm=[]; + for ki=1:p + Rm=[Rm; Rk(:,:,ki)]; + end + R(:,:,k)=A*Rm; + end + +end + diff --git a/toolbox/connectivity/YuleWalker_Mask.m b/toolbox/connectivity/YuleWalker_Mask.m new file mode 100644 index 000000000..eca8525e7 --- /dev/null +++ b/toolbox/connectivity/YuleWalker_Mask.m @@ -0,0 +1,38 @@ +function [A,Sigma] = YuleWalker_Mask(Gamma, A_mask) + +[N,~,q] = size(Gamma); +p = q - 1; + +% Creo le matrici che mi servono per il conto +Psi = zeros(N*p); +G = reshape(Gamma(:,:,2:end),[N,N*p]); +% B_mask = reshape(B_mask,[N,N*p]); +A_mask = repmat(A_mask,[1,p]); + +for i = 1:p + for j = 1:p + k = j - i; + if (k > 0) % Colonna > Riga = metà superiore destra + Psi((i-1)*N+1:i*N, (j-1)*N+1:j*N) = Gamma(:,:,k + 1); + else + Psi((i-1)*N+1:i*N, (j-1)*N+1:j*N) = Gamma(:,:,abs(k) + 1)'; + end + end +end + +% Mi muovo riga per riga +A = zeros(N,N*p); + +for i = 1:N + % Vedo nella maschera quali elementi fisso a zero + ind = find(A_mask(i,:) ~= 0); + Psi_t = Psi(ind,ind); + G_t = G(i,ind); + + A(i,ind) = G_t / Psi_t; +end + +Sigma = Gamma(:,:,1) - A * G'; + +end + diff --git a/toolbox/connectivity/bst_granger.m b/toolbox/connectivity/bst_granger.m index e3fa4b4d4..d8393798e 100644 --- a/toolbox/connectivity/bst_granger.m +++ b/toolbox/connectivity/bst_granger.m @@ -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 @@ -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: @@ -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) @@ -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 @@ -153,17 +116,10 @@ %% 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 -end - if isempty(Y) % auto-causality % setup connectivity = zeros(nX); - connectivityV = zeros(nX); % only iterate over one triangle for iX = 1:nX @@ -171,57 +127,37 @@ % 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 - - % 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 = YuleWalker_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] = YuleWalker_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 - - % 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 - catch ME %#ok % 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] = YuleWalker_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 @@ -233,7 +169,6 @@ % setup connectivity = zeros(nX, nY); - connectivityV = zeros(nX, nY); duplicates = zeros(0, 2); for iX = 1:nX @@ -241,45 +176,21 @@ 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 - - % 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 - - % 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 - catch ME %#ok % 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 = YuleWalker_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] = YuleWalker_Mask(R, mask); + + % connectivity + connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); else % save duplicates to modify later @@ -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) + \ No newline at end of file diff --git a/toolbox/connectivity/bst_granger_conditional.m b/toolbox/connectivity/bst_granger_conditional.m new file mode 100644 index 000000000..25ffc8e83 --- /dev/null +++ b/toolbox/connectivity/bst_granger_conditional.m @@ -0,0 +1,236 @@ +function [connectivity, pValue] = bst_granger_conditional(X, Y, order, inputs) +% BST_GRANGER_CONDITIONAL Conditional Granger causality between any +% two signals +% +% Inputs: +% sinks - first set of signals, one signal per row +% [X: MX x N or MX x N x T matrix] +% sources - second set of signals, one signal per row +% [Y: MY x N or MY x N x T matrix] +% (default: Y = X) +% order - maximum lag in AR model for causality in mean +% [p: nonnegative integer] +% inputs - structure of parameters: +% |-nTrials - # of trials in concantenated signal +% | [T: positive integer] +% |-standardize - if true (default), remove mean from each signal. +% | if false, assume signal has already been detrended +% |-flagFPE - if true, optimize order for AR model +% | if false (default), force same order in all AR models +% | [E: default false] +% +% Outputs: +% connectivity - A x B matrix of causalities from source to sink +% [C: MX x MY matrix] +% pValue - parametric p-value for corresponding conditional Granger causality estimate +% [P: MX x MY matrix] (not implemented!) +% +% See also BST_MVAR +% +% For each signal pair (a,b) we calculate the conditional Granger causality: +% +% det(restricted_residuals_cov_matrix) +% gc = ---------------------------------------------------- +% det(full_cov_matrix) +% +% see Barret, Barnett and Seth "Multivariate Granger causality and generalized variance". +% +% Call: +% connectivity = bst_granger_conditional(X, Y, 5, inputs); +% connectivity = bst_granger_conditional(X, [], 5, inputs); % every pair in X +% connectivity = bst_granger_conditional(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 + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Sergul Aydore & Syed Ashrafulla, 2012 +% Modified by: Davide Nuzzi, 2021 + +% default: 1 trial +if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) + inputs.nTrials = 1; +end + +% default: do not optimize order in MVAR modeling +if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) + inputs.flagFPE = false; +end + +% dimensions of the signals +nX = size(X, 1); +if ndims(X) == 3 + inputs.nTrials = size(X,3); + X = reshape(X, nX, []); +end +nSamples = size(X, 2); +nTimes = nSamples / inputs.nTrials; + +% if we are doing auto-causality, empty Y so that we only look at X +if ~isempty(Y) + nY = size(Y, 1); + Y = reshape(Y, nY, []); +end + +% remove linear effects if desired +if isfield(inputs, 'standardize') && inputs.standardize + detrender = [ ... + ones(1, nTimes); ... % constant trend in data + linspace(-1, 1, nTimes); ... % linear trend in data + 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data + ]; + + % detrend X + for iTrial = 1:inputs.nTrials + X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); + end + + % detrend Y only if it is not empty + if ~isempty(Y) + for iTrial = 1:inputs.nTrials + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); + end + end +end + +%% Iterate over all pairs of sinks & sources + +if isempty(Y) % auto-causality + + % setup + connectivity = zeros(nX); + + % the full model is evaluated only one time + + % multivariate model estimation + [transfers, noiseCovariance, order] = bst_mvar(X, order, inputs.nTrials, inputs.flagFPE); + + % data correlations using Yule-Walker (up to high order 50) + R = YuleWalker_Inverse(transfers, noiseCovariance, 50); + + % only iterate over one triangle + for iX = 1:nX + + % iterate over all the pairs after iX + for iY = (iX+1):nX + + % restricted model iY -> iX + + % mask for the coefficients of the restricted model + mask = ones(size(X,1)); + mask(iX,iY) = 0; + + % restricted multivariate model (using masked row-by-row solution of YW equations) + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); + + % connectivity + connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); + + % restricted model iX -> iY + + % mask for the coefficients of the restricted model + mask = ones(size(X,1)); + mask(iY,iX) = 0; + + % restricted bivariate model (using masked row-by-row solution of + % YW equations) + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_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 + connectivity(iX, iX) = max([connectivity(iX, :) connectivity(:, iX)']); + + end + +else % cross-causality + + % setup + connectivity = zeros(nX, nY); + duplicates = zeros(0, 2); + + for iX = 1:nX + for iY = 1:nY + % X(iX,:) = sink, Y(iY,:) = source, Y(~iY,:) = conditional + % I have to check that the sink is different form all the other + % variables in Y + sink = X(iX,:); + source = Y(iY,:); + other_vars = Y(setdiff( 1:size(Y, 1), iY),:); + + % If sink and source are duplicates there's a correction + if max(abs(sink - source)) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere + + % If the target is equal to one of the conditioning + % variables, remove it + remove_inds = []; + for k = 1:size(other_vars,1) + if max(abs(sink - other_vars(k,:))) < eps + remove_inds = [remove_inds, k]; + end + end + other_vars(remove_inds,:) = []; + + + % multivariate model estimation (one sink, all sources) + [transfers, noiseCovariance, order] = bst_mvar([sink; source; other_vars], order, inputs.nTrials, inputs.flagFPE); + + % data correlations using Yule-Walker (up to high order 50) + R = YuleWalker_Inverse(transfers, noiseCovariance, 50); + + % mask for the coefficients of the restricted model + mask = ones(2 + size(other_vars,1)); % size(other vars) + one source + one sink + mask(1,2) = 0; % Cut only the source -> sink coupling + + % restricted multivariate model (using masked row-by-row solution of YW equations) + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); + + % connectivity + connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); + + else % save duplicates to modify later + + duplicates(end+1, :) = [iX iY]; %#ok + + end + + end + end + + % for duplicate indices, set the causality value to the maximum over all inflows for iX and outflows for iY + for iDuplicate = 1:size(duplicates, 1) + connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2)) = ... + max([connectivity(duplicates(iDuplicate, 1), :) connectivity(:, duplicates(iDuplicate, 2))']); + end + +end + +%% Statistics: parametric p-values for causality in mean (based on regression coefficients) and variance (based on Wald statistics) + +% F statistic of connectivity when multiplied by number of regressors +pValue = 1; +% 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) + \ No newline at end of file diff --git a/toolbox/connectivity/bst_granger_spectral.m b/toolbox/connectivity/bst_granger_spectral.m index 6f2168894..b07ed5c3f 100644 --- a/toolbox/connectivity/bst_granger_spectral.m +++ b/toolbox/connectivity/bst_granger_spectral.m @@ -24,44 +24,35 @@ % % Outputs: % connectivity - A x B matrix of spectral Granger causalities from -% source to sink. For each signal pair (a,b) we calculate -% -% S_{sink} (f) -% ------------------------------------------------------------ -% S_{sink}(f) - |H_{sink, source} (f)|^2 sigma_{source | sink} -% -% with S_{sink}(f) as the power spectral density of a @ f -% H_{sink, source}(f) as the transfer function @ f -% sigma_{source | sink} as the conditional variance -% of the residual at b given the residual at a, -% calculated using the residual covariance. +% source to sink. % By default, GC(a,a) = 0 if Y is empty. % [C: MX x MY x NF matrix] % pValues - parametric p-value for corresponding spectral Granger -% causality in mean estimate +% causality in mean estimate (not implemented!) % [P: MX x MY x NF matrix] % freq - frequencies corresponding to the previous two metrics % [F: NF x 1 vector] % -% See also BST_GRANGER, BST_COHERENCE_MVAR. +% See also BST_GRANGER. +% +% Spectral causality measures are evaluated from the spectra of the full +% and restricted models as: +% +% gc(w) = det(S_restricted(w))/det(S_full(w)) +% +% see Cohen, Dror, et al. "A general spectral decomposition of causal influences +% applied to integrated information." Journal of neuroscience methods 330 (2020): 108443 +% for additional information. % % Call: % connectivity = bst_granger_spectral(X, Y, 200, 10, inputs); % general call % connectivity = bst_granger_spectral(X, [], 200, 30, inputs); % every pair -% connectivity = bst_granger_spectral(X, [], 200, 30, inputs); % more variance % Parameter examples: % inputs.freq = 0:0.1:100; % specify desired frequencies % inputs.freqResolution = 0.1; % have a high-point FFT % inputs.nTrials = 9; % use trial-average covariances in AR estimation % inputs.flagFPE = true; % use AR model with best information criteria -% Note: for those following Chicharro2012, I have used the equivalence -% sigma_{xx}^(xy) |H_{xx}^(xy) (w)|^2 -% = -% S_{xx}(w) - sigma_{yy}^(xy) |H_{xy}^(xy) (w)|^2 -% which follows Geweke1982 instead. As Chicharro notes, it may be better to use his -% formulation because it separates out instantaneous causality I think. - % @============================================================================= % This function is part of the Brainstorm software: % https://neuroimage.usc.edu/brainstorm @@ -81,11 +72,13 @@ % =============================================================================@ % % Authors: Sergul Aydore & Syed Ashrafulla, 2012 +% Modified by: Davide Nuzzi, 2021 % reformat to a 2D matrix if necessary, and pull out # of trials if ndims(X) == 3 inputs.nTrials = size(X,3); X = reshape(X, size(X, 1), []); + Y = reshape(Y, size(Y, 1), []); elseif ~isfield(inputs, 'nTrials') inputs.nTrials = 1; end @@ -110,7 +103,6 @@ % detrend Y only if it is not empty if ~isempty(Y) - Y = reshape(Y, size(Y, 1), []); % reshape to 2D matrix first for iTrial = 1:inputs.nTrials Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); @@ -157,26 +149,49 @@ for iX = 1:size(X, 1) for iY = iX+1 : size(X, 1) % to avoid auto-causality - % two-variate model for given source + % two-variate model [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); X(iY, :)], order, inputs.nTrials, inputs.flagFPE); + + % spectra of full model + [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); - % spectra and power of forward system - [spectra, freq, forward] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); + % data correlations using Yule-Walker (up to high order 50) + R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - % Geweke-Granger spectral causality from source to sink - unrestricted = abs(spectra(1, 1, :)); restriction = forward(1, 2, :); % S_{sink}(f) and |H_{sink, source} (f)|^2, w/ abs to get rid of 1e-16 imag part - residualVariance = noiseCovariance(2,2) - noiseCovariance(2, 1) / noiseCovariance(1, 1) * noiseCovariance(1, 2); % partial variance of source - restricted = abs(unrestricted) - abs(restriction).^2 * residualVariance; % S_{sink} (f) - |H_{sink, source} (f)|^2 sigma_{source | sink} - connectivity(iX, iY, abs(restricted) > 1e-60) = unrestricted(abs(restricted) > 1e-60) ./ restricted(abs(restricted) > 1e-60) - 1;% Geweke-Granger - % sigma_{source | sink} = partial covariance which is the formula above (sigma_{source} - rho_{source, sink} / sigma_{sink} * rho_{sink, source}) + % restricted model iY -> iX + + % 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] = YuleWalker_Mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + + % restricted model iX -> iY - % Geweke-Granger spectral causality from sink back to source (to halve the # of MVAR fittings) - unrestricted = abs(spectra(2, 2, :)); restriction = forward(2, 1, :); % S_{source}(f) and |H_{source, sink} (f)|^2, w/ abs to get rid of 1e-16 imag part - residualVariance = noiseCovariance(1,1) - noiseCovariance(1, 2) / noiseCovariance(2, 2) * noiseCovariance(2, 1); % partial variance of sink - restricted = abs(unrestricted) - abs(restriction).^2 * residualVariance; % S_{source} (f) - |H_{source, sink} (f)|^2 sigma_{sink | source} - connectivity(iY, iX, abs(restricted) > 1e-60) = unrestricted(abs(restricted) > 1e-60) ./ restricted(abs(restricted) > 1e-60) - 1; % Geweke-Granger - % sigma_{sink | source} = partial covariance which is the formula above (sigma_{sink} - rho_{sink, source} / sigma_{source} * rho_{source, sink}) + % mask for the coefficients of the restricted model + mask = ones(2); + mask(2,1) = 0; + % restricted bivariate model + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iY, iX, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end end % diagonal will equal the maximum of all inflows and outflows for iX, specific to each frequency @@ -196,22 +211,30 @@ if max(abs(X(iX, :) - Y(iY, :))) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere - % 2-variate model for given source + % two-variate model [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); Y(iY, :)], order, inputs.nTrials, inputs.flagFPE); - - % spectra and power of forward system - [spectra, freq, forward] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); - spectra = abs(spectra(1, 1, :)); % limit to the autospectrum of the sink S_{sink} (f) and take absolute value to rid the 1e-16 imaginary part - forward = forward(1, 2, :); % limit to the cross-transfer from source to sink H_{sink, source} (f) - % partial covariance of residual - residualVariance = noiseCovariance(2,2) - noiseCovariance(2,1) / noiseCovariance(1, 1) * noiseCovariance(1, 2); - % sigma_{source | sink} = partial covariance which is the formula above (sigma_{source} - rho_{source, sink} / sigma_{sink} * rho_{sink, source}) + % spectra of full model + [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); + + % data correlations using Yule-Walker (up to high order 50) + R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - % Geweke-Granger spectral causality from source to sink - restricted = spectra - abs(forward).^2 * residualVariance; % S_{sink} (f) - |H_{sink, source} (f)|^2 sigma_{source | sink} - connectivity(iX, iY, abs(restricted) > 1e-60) = spectra(abs(restricted) > 1e-60) ./ restricted(abs(restricted) > 1e-60) - 1; % Geweke-Granger + % 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] = YuleWalker_Mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + else % save duplicates to modify later duplicates(end+1, :) = [iX iY]; %#ok @@ -243,18 +266,18 @@ end -%% ======================================================== estimation for multivariate autoregression ======================================================== -function [spectra, freq, forward] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs) -% BST_MVAR_SPECTRUM Calculate the parametric spectra of a bivariate system -% with given MVAR coefficients & an estimate of the -% covariance matrix of the innovation (i.e. noise) -% process. +%% ======================================================== spectra estimation ======================================================== +function [spectra, freq] = bst_granger_spectral_spectrum(A, Sigma, nFFT, Fs) +% BST_GRANGER_SPECTRAL_SPECTRUM Calculate the parametric spectra of a multivariate system +% with given MVAR coefficients and an estimate of the +% covariance matrix of the innovation (i.e. noise) +% process. % % Inputs: % transfers - transfer matrices in AR process -% [A: 2 x 2P matrix, P = order] +% [A: N x NP matrix, P = order] % noiseCovariance - variance of residuals -% [C: 2 x 2 matrix] +% [C: N x N matrix] % nFFT - number of FFT bins to calculate spectra % [NF: positive number, usually power of 2] % Fs - sampling frequency of data @@ -262,84 +285,40 @@ % % Outputs: % spectra - cross-spectrum between each pair of variables -% [S: 2 x 2 x NF/2 matrix] +% [S: N x N x NF/2 matrix] % freq - frequencies used based on # of FFT bins % [F: NF/2 x 1 matrix] -% forward - forward transform in frequency from source to sink -% [H: 2 x 2 x NF/2 matrix] % --> all outputs have length NF/2, ignoring frequenices past Fs/2 -% % Call: % spectra = bst_mvar_spectrum(transfers, C, 512, 200); % basic % spectra = bst_mvar_spectrum(transfers, C, 2048, 200); % add FFT interp -% [spectra, freq] = bst_mvar_spectrum(transfers, C, 512, 200); % grab freqs -% [~, ~, forward] = bst_mvar_spectrum(transfers, C, 64, 200); % causality - -% Notes: -% The DTFT we want is -% H = I - sum_p C_p e^{-j 2 pi f/Fs p} = sum_p D_p e^{-j 2 pi f/Fs (p-1)} where D_1 = 1 and D_p = -C_{p-1} for p > 1 -% MatLab's FFT provides -% G = sum_n x_n e^{-j 2 pi (k-1)/N (n-1)} -% so the matching is p = n and (k-1)/N = f/Fs, after vectorizing H. % frequencies to estimate cross-spectra freq = Fs/2*linspace(0, 1, nFFT/2 + 1); freq(end) = []; -%% Inverse of transfer function in frequency -% Inverse transfer means the transfer function from the sources to the innovations. We calculate this at each frequency. -% The inverse transfer from source a to innovation b is -% 1 - \sum_{n=1}^N a_{ab} [n] e^{-j2pi * f * n} - -inverse = fft(reshape([eye(2) -transfers], 4, [])', nFFT); % reshape so we can do a vector FFT quickly -inverse = inverse(1:nFFT/2, :); % restrict to the first symmetric half of the spectrum +% get number of variables, fft bins and order of the model +M = nFFT/2; +[N,pN] = size(A); +p = pN/N; -%% Forward transfer in autoregressive model +% evaluate the trasfer function and the spectra +H = complex(zeros(N,N,M)); +spectra = complex(zeros(N,N,M)); -% pieces of 2x2 inverse -forward = zeros(2, 2, nFFT/2); % an important caveat is that I did not reshape inverse earlier; -forward(1,1,:) = inverse(:,4); forward(1,2,:) = -inverse(:,3); % as a result, we have a column-wise index mapping: 4 = (2,2) and 3 = (1,2) -forward(2,1,:) = -inverse(:,2); forward(2,2,:) = inverse(:,1); % in addition, 2 = (2,1) and 1 = (1,1). then these elements fit the 2x2 matrix inverse -detInverse = inverse(:,1).*inverse(:,4) - inverse(:,3).*inverse(:,2); % the same thing happens here, using the indexing to avoid a reshape() call +A_r = reshape(A,[N,N,p]); +w_vec = 0:pi/(nFFT/2-1):pi; -% normalization by determinant to get inverse -forward = bst_bsxfun(@rdivide, forward, reshape(detInverse, [1 1 length(freq)])); % complete the inversion by dividing by frequency-dependent determinant - -%% Cross-spectrum from forward model -% The forward transfer is the inverse of the inverse transfer at each frequency f. Denoted H(f), the power spectral density is then HH' at each frequency. - -% the loop that we won't use -% for idxFreq = 1:nFFT -% spectra(:, :, idxFreq) = H(:,:,idxFreq) * noiseCovariance * H(:,:,idxFreq)'; -% end - -% formula for the matrix multiplication -spectra(1,1,:) = ... % 1,1 element is H_11 C_11 H_11^* + H_12 C_12 H_11^* + H_11 C_12 H_12^* + H_12 C_22 H_12^* (and I combine the middle two into 2 Re{.}) - noiseCovariance(1,1) * forward(1,1,:) .* conj(forward(1,1,:)) ... - + noiseCovariance(1,2) * real(forward(1,2,:) .* conj(forward(1,1,:))) * 2 ... - + noiseCovariance(2,2) * forward(1,2,:) .* conj(forward(1,2,:)); -spectra(1,2,:) = ... % 1,2 element is H_11 C_11 H_21^* + H_12 C_12 H_21^* + H_11 C_12 H_22^* + H_12 C_22 H_22^* - noiseCovariance(1,1) * forward(2,1,:) .* conj(forward(1,1,:)) ... - + noiseCovariance(1,2) * forward(2,2,:) .* conj(forward(1,1,:)) ... - + noiseCovariance(1,2) * forward(2,1,:) .* conj(forward(1,2,:)) ... - + noiseCovariance(2,2) * forward(1,2,:) .* conj(forward(1,2,:)); -spectra(2,1,:) = conj(spectra(1,2,:)); % for speed, force the 2,1 element to be the conjugate of the 1,2 element so we have conjugate symmetry -spectra(2,2,:) = ... % 2,2 element is H_21 C_11 H_21^* + H_22 C_12 H_21^* + H_21 C_12 H_22^* + H_22 C_22 H_22^* (and I combine the middle two into 2 Re{.}) - noiseCovariance(1,1) * forward(2,1,:) .* conj(forward(2,1,:)) ... - + noiseCovariance(1,2) * real(forward(2,2,:) .* conj(forward(2,1,:))) * 2 ... - + noiseCovariance(2,2) * forward(2,2,:) .* conj(forward(2,2,:)); - -%% Normalize for sampling rate, and truncate if necessary - -% normalization -forward = forward / sqrt(Fs); -spectra = spectra / Fs; +for n = 1:M + e = zeros(p,1); + for k = 1:p + e(k) = exp(-1i * w_vec(n) * k); + end + e = permute(repmat(e,[1,N,N]),[2,3,1]); + A_w = eye(N) - sum(A_r .* e,3); + H(:,:,n) = inv(A_w); + spectra(:,:,n) = H(:,:,n) * Sigma * ctranspose(H(:,:,n)); +end -%% Confidence intervals -% taken from Kay, pp 194-195 -% alpha = 1 - ci; -% original = -sqrt(2) * erfcinv( 2 * ( 1 - alpha/2) ); -% lower = abs(spectra) * (1 - sqrt(2 * order / nTimes) * original); -% upper = abs(spectra) * (1 + sqrt(2 * order / nTimes) * original); end %% <== FUNCTION END \ No newline at end of file diff --git a/toolbox/connectivity/bst_granger_spectral_conditional.m b/toolbox/connectivity/bst_granger_spectral_conditional.m new file mode 100644 index 000000000..f2015044b --- /dev/null +++ b/toolbox/connectivity/bst_granger_spectral_conditional.m @@ -0,0 +1,350 @@ +function [connectivity, pValues, freq] = bst_granger_spectral_conditional(X, Y, Fs, order, inputs) +% BST_GRANGER_SPECTRAL_CONDITIONAL Granger causality at each frequency between any two +% signals, conditioned on all the remaining signals +% +% Inputs: +% X - first set of signals, one signal per row +% [X: A x N or A x N x T matrix] +% Y - second set of signals, one signal per row +% [Y: B x N or B x N x T matrix] +% (default: Y = X) +% Fs - sampling rate (we assume uniform sampling rate) +% [FS: scalar, FS > freq(end)*2] +% order - maximum order of autogressive model +% [p: integer > 1, default 10] +% inputs - structure of parameters: +% |-freq - frequencies of interest if desired +% |-freqResolution - maximum freq resolution in Hz, to limit NFFT +% | [DF: double, default [] (i.e. no limit)] +% |-nTrials - # of trials in concantenated signal +% |-flagFPE - if true, optimize order for autoregression +% | if false (default), use same order in autoregression +% |-standardize - if true (default), remove mean from each signal +% | if false, assume signal has already been detrended +% +% Outputs: +% connectivity - A x B matrix of spectral Granger causalities from +% source to sink, conditioned on all the other variables. +% If Y is empty, the matrix contains the spectral GC +% from each source variable to each sink variable, +% conditioned on the remaining variable. +% If Y is not empty the matrix contains the spectral +% causalities from each variable in Y to each +% variable in X, conditionend on the other variables +% in Y. +% [C: MX x MY x NF matrix] +% pValues - parametric p-value for corresponding spectral Granger +% causality in mean estimate (not implemented yet!) +% [P: MX x MY x NF matrix] +% freq - frequencies corresponding to the previous two metrics +% [F: NF x 1 vector] +% +% See also BST_GRANGER_CONDITIONAL. +% + +% Spectral causality measures are evaluated from the spectra of the full +% and restricted models as: +% +% gc(w) = det(S_restricted(w))/det(S_full(w)) +% +% see Cohen, Dror, et al. "A general spectral decomposition of causal influences +% applied to integrated information." Journal of neuroscience methods 330 (2020): 108443 +% for additional information. +% +% Call: +% connectivity = bst_granger_spectral_conditional(X, Y, 200, 10, inputs); % general call +% connectivity = bst_granger_spectral_conditional(X, [], 200, 30, inputs); % every pair in X +% Parameter examples: +% inputs.freq = 0:0.1:100; % specify desired frequencies +% inputs.freqResolution = 0.1; % have a high-point FFT +% inputs.nTrials = 9; % use trial-average covariances in AR estimation +% inputs.flagFPE = true; % use AR model with best information criteria + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Sergul Aydore & Syed Ashrafulla, 2012 +% Modified by: Davide Nuzzi, 2021 + +% reformat to a 2D matrix if necessary, and pull out # of trials +if ndims(X) == 3 + inputs.nTrials = size(X,3); + X = reshape(X, size(X, 1), []); + Y = reshape(Y, size(Y, 1), []); +elseif ~isfield(inputs, 'nTrials') + inputs.nTrials = 1; +end + +% lengths of things +nSamples = size(X, 2); +nTimes = nSamples / inputs.nTrials; + +% standardization: zero mean & unit variance, remove linear & quadratic trends as well +if isfield(inputs, 'standardize') && inputs.standardize + detrender = [ ... + ones(1, nTimes); ... % constant trend in data + linspace(-1, 1, nTimes); ... % linear trend in data + 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data + ]; + + % detrend X + for iTrial = 1:inputs.nTrials + X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); + end + + % detrend Y only if it is not empty + if ~isempty(Y) + for iTrial = 1:inputs.nTrials + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); + end + end +end + +% number of FFT bins required +if ~isfield(inputs, 'freq') || isempty(inputs.freq) % frequencies of interest are not defined + if isfield(inputs, 'freqResolution') && ~isempty(inputs.freqResolution) && (size(X,2) > round(Fs / inputs.freqResolution)) + nFFT = 2^nextpow2( round(Fs / inputs.freqResolution) ); + else % use a default frequency resolution of 1Hz to mirror the standard frequency resolution in bst_coherence_welch.m + nFFT = 2^nextpow2( round(Fs / 1) ); + end +elseif numel(inputs.freq) == 1 + nFFT = inputs.freq; +else + nFFT = 2^nextpow2(max( length(inputs.freq)-1, (Fs/2) / min(diff( sort(inputs.freq(:)) )) )) * 2; +end + +% default: Order 10 used in BrainStorm +if isempty(order) + order = 10; +end + +% default: Single-trial +if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) + inputs.nTrials = 1; +end + +% default: do not optimize model order +if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) + inputs.flagFPE = false; +end + +%% Differentiate between auto-causality and cross-causality for speed + +if isempty(Y) % causality between signals in X + + % preallocate the spectral causality matrix + connectivity = zeros(size(X, 1), size(X, 1), nFFT/2); + + % the full model is evaluated only one time + + % multivariate model estimation + [transfers, noiseCovariance, order] = bst_mvar(X, order, inputs.nTrials, inputs.flagFPE); + + % spectra for the full model + [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); + + % data correlations using Yule-Walker (up to high order 50) + R = YuleWalker_Inverse(transfers, noiseCovariance, 50); + + % iterate over all pairs of sinks & sources + for iX = 1:size(X, 1) + for iY = iX+1 : size(X, 1) % to avoid auto-causality + + % restricted model iY -> iX + + % mask for the coefficients of the restricted model + mask = ones(size(X,1)); + mask(iX,iY) = 0; + + % restricted multivariate model (using masked row-by-row solution of YW equations) + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + + % restricted model iX -> iY + + % mask for the coefficients of the restricted model + mask = ones(size(X,1)); + mask(iY,iX) = 0; + + % restricted multivariate model (using masked row-by-row solution of YW equations) + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iY, iX, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + end + + % diagonal will equal the maximum of all inflows and outflows for iX, specific to each frequency + connectivity(iX, iX, :) = max( max(connectivity(iX, :, :), [], 2), max(connectivity(:, iX, :), [], 1) ); + end + +else % we have to use all pairs of signals + + % preallocate the spectral causality matrix + connectivity = zeros(size(X, 1), size(Y, 1), nFFT/2); + duplicates = zeros(0, 2); + + % iterate over all pairs of sinks & sources + for iX = 1:size(X, 1) + for iY = 1:size(Y, 1) + + % X(iX,:) = sink, Y(iY,:) = source, Y(~iY,:) = conditional + % I have to check that the sink is different form all the other + % variables in Y + sink = X(iX,:); + source = Y(iY,:); + other_vars = Y(setdiff( 1:size(Y, 1), iY),:); + + % If sink and source are duplicates there's a correction + % + if max(abs(sink - source)) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere + + % If the target is equal to one of the conditioning + % variables, remove it + remove_inds = []; + for k = 1:size(other_vars,1) + if max(abs(sink - other_vars(k,:))) < eps + remove_inds = [remove_inds, k]; + end + end + other_vars(remove_inds,:) = []; + + + % multivariate model estimation (one sink, all sources) + [transfers, noiseCovariance, order] = bst_mvar([sink; source; other_vars], order, inputs.nTrials, inputs.flagFPE); + + % spectra for the full model + [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); + + % data correlations using Yule-Walker (up to high order 50) + R = YuleWalker_Inverse(transfers, noiseCovariance, 50); + + % mask for the coefficients of the restricted model + mask = ones(2 + size(other_vars,1)); % size(other vars) + one source + one sink + mask(1,2) = 0; % Cut only the source -> sink coupling + + % restricted multivariate model (using masked row-by-row solution of YW equations) + [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + + else % save duplicates to modify later + + duplicates(end+1, :) = [iX iY]; %#ok + + end + + end + end + + % for duplicate indices, set the causality value to the maximum of all inflows for iX and outflows for iY + for iDuplicate = 1:size(duplicates, 1) + connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2), :) = max( ... + max(connectivity(duplicates(iDuplicate, 1), :, :), [], 2), ... + max(connectivity(:, duplicates(iDuplicate, 2), :), [], 1) ... + ); + end + +end + +%% Interpolate to desired frequencies & perform statistics if desired +if isfield(inputs, 'freq') && ~isempty(inputs.freq) && numel(inputs.freq) > 1 + connectivity = permute(interp1(freq, permute(connectivity, [3 1 2]), inputs.freq), [2 3 1]); + % pValues = permute(interp1(freq, permute(pValues, [3 1 2]), inputs.freq), [2 3 1]); + freq = inputs.freq; +end + +pValues = NaN; % no parametric p-values for now + +end + +%% ======================================================== spectra estimation ======================================================== +function [spectra, freq] = bst_granger_spectral_spectrum(A, Sigma, nFFT, Fs) +% BST_GRANGER_SPECTRAL_SPECTRUM Calculate the parametric spectra of a multivariate system +% with given MVAR coefficients and an estimate of the +% covariance matrix of the innovation (i.e. noise) +% process. +% +% Inputs: +% transfers - transfer matrices in AR process +% [A: N x NP matrix, P = order] +% noiseCovariance - variance of residuals +% [C: N x N matrix] +% nFFT - number of FFT bins to calculate spectra +% [NF: positive number, usually power of 2] +% Fs - sampling frequency of data +% [FS: double, FS > freq(end)*2] +% +% Outputs: +% spectra - cross-spectrum between each pair of variables +% [S: N x N x NF/2 matrix] +% freq - frequencies used based on # of FFT bins +% [F: NF/2 x 1 matrix] +% --> all outputs have length NF/2, ignoring frequenices past Fs/2 +% Call: +% spectra = bst_mvar_spectrum(transfers, C, 512, 200); % basic +% spectra = bst_mvar_spectrum(transfers, C, 2048, 200); % add FFT interp + +% frequencies to estimate cross-spectra +freq = Fs/2*linspace(0, 1, nFFT/2 + 1); +freq(end) = []; + +% get number of variables, fft bins and order of the model +M = nFFT/2; +[N,pN] = size(A); +p = pN/N; + +% evaluate the trasfer function and the spectra +H = complex(zeros(N,N,M)); +spectra = complex(zeros(N,N,M)); + +A_r = reshape(A,[N,N,p]); +w_vec = 0:pi/(nFFT/2-1):pi; + +for n = 1:M + e = zeros(p,1); + for k = 1:p + e(k) = exp(-1i * w_vec(n) * k); + end + e = permute(repmat(e,[1,N,N]),[2,3,1]); + A_w = eye(N) - sum(A_r .* e,3); + H(:,:,n) = inv(A_w); + spectra(:,:,n) = H(:,:,n) * Sigma * ctranspose(H(:,:,n)); +end + + +end %% <== FUNCTION END \ No newline at end of file From b5d3f1ecc029a88294d316cf0d4c122c4b1ce72c Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 18:58:42 +0200 Subject: [PATCH 2/9] Delete YuleWalker_Mask.m --- toolbox/connectivity/YuleWalker_Mask.m | 38 -------------------------- 1 file changed, 38 deletions(-) delete mode 100644 toolbox/connectivity/YuleWalker_Mask.m diff --git a/toolbox/connectivity/YuleWalker_Mask.m b/toolbox/connectivity/YuleWalker_Mask.m deleted file mode 100644 index eca8525e7..000000000 --- a/toolbox/connectivity/YuleWalker_Mask.m +++ /dev/null @@ -1,38 +0,0 @@ -function [A,Sigma] = YuleWalker_Mask(Gamma, A_mask) - -[N,~,q] = size(Gamma); -p = q - 1; - -% Creo le matrici che mi servono per il conto -Psi = zeros(N*p); -G = reshape(Gamma(:,:,2:end),[N,N*p]); -% B_mask = reshape(B_mask,[N,N*p]); -A_mask = repmat(A_mask,[1,p]); - -for i = 1:p - for j = 1:p - k = j - i; - if (k > 0) % Colonna > Riga = metà superiore destra - Psi((i-1)*N+1:i*N, (j-1)*N+1:j*N) = Gamma(:,:,k + 1); - else - Psi((i-1)*N+1:i*N, (j-1)*N+1:j*N) = Gamma(:,:,abs(k) + 1)'; - end - end -end - -% Mi muovo riga per riga -A = zeros(N,N*p); - -for i = 1:N - % Vedo nella maschera quali elementi fisso a zero - ind = find(A_mask(i,:) ~= 0); - Psi_t = Psi(ind,ind); - G_t = G(i,ind); - - A(i,ind) = G_t / Psi_t; -end - -Sigma = Gamma(:,:,1) - A * G'; - -end - From 131c1e2b486fbfb03648821b0472ab5dd9c04c65 Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 18:58:48 +0200 Subject: [PATCH 3/9] Delete YuleWalker_Inverse.m --- toolbox/connectivity/YuleWalker_Inverse.m | 34 ----------------------- 1 file changed, 34 deletions(-) delete mode 100644 toolbox/connectivity/YuleWalker_Inverse.m diff --git a/toolbox/connectivity/YuleWalker_Inverse.m b/toolbox/connectivity/YuleWalker_Inverse.m deleted file mode 100644 index 48def697f..000000000 --- a/toolbox/connectivity/YuleWalker_Inverse.m +++ /dev/null @@ -1,34 +0,0 @@ -function R = YuleWalker_Inverse(A,Sigma, maxOrder) -%YULEWALKER Summary of this function goes here -% A partire dalla matrice A e dalle correlazioni dei residui ottengo le -% correlazioni a vari ritardi. Le Gamma sono matrici N x N e in totale -% saranno (P+1) -% L'approccio dovrebbe essere di trasformare il problema in un'equazione -% di Lyapunov - - [N,~] = size(Sigma); - p = size(A,2)/N; - - % Copiato da CElinVAR - Im = eye(N*p); - F = [A; Im(1:end-N,:)]; - Delta=zeros(N*p); - Delta(1:N,1:N) = Sigma; - BigSigma = dlyap(F,Delta); - - R=NaN*ones(N,N,maxOrder+1); - for i=1:p - R(:,:,i)=BigSigma(1:N,N*(i-1)+1:N*i); - end - - for k=p+1:maxOrder+1 - Rk=R(:,:,k-1:-1:k-p); - Rm=[]; - for ki=1:p - Rm=[Rm; Rk(:,:,ki)]; - end - R(:,:,k)=A*Rm; - end - -end - From 17b3c57c023825da5d7729cb4d4b4fbe4525e6fe Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 18:59:03 +0200 Subject: [PATCH 4/9] Delete bst_granger_conditional.m --- .../connectivity/bst_granger_conditional.m | 236 ------------------ 1 file changed, 236 deletions(-) delete mode 100644 toolbox/connectivity/bst_granger_conditional.m diff --git a/toolbox/connectivity/bst_granger_conditional.m b/toolbox/connectivity/bst_granger_conditional.m deleted file mode 100644 index 25ffc8e83..000000000 --- a/toolbox/connectivity/bst_granger_conditional.m +++ /dev/null @@ -1,236 +0,0 @@ -function [connectivity, pValue] = bst_granger_conditional(X, Y, order, inputs) -% BST_GRANGER_CONDITIONAL Conditional Granger causality between any -% two signals -% -% Inputs: -% sinks - first set of signals, one signal per row -% [X: MX x N or MX x N x T matrix] -% sources - second set of signals, one signal per row -% [Y: MY x N or MY x N x T matrix] -% (default: Y = X) -% order - maximum lag in AR model for causality in mean -% [p: nonnegative integer] -% inputs - structure of parameters: -% |-nTrials - # of trials in concantenated signal -% | [T: positive integer] -% |-standardize - if true (default), remove mean from each signal. -% | if false, assume signal has already been detrended -% |-flagFPE - if true, optimize order for AR model -% | if false (default), force same order in all AR models -% | [E: default false] -% -% Outputs: -% connectivity - A x B matrix of causalities from source to sink -% [C: MX x MY matrix] -% pValue - parametric p-value for corresponding conditional Granger causality estimate -% [P: MX x MY matrix] (not implemented!) -% -% See also BST_MVAR -% -% For each signal pair (a,b) we calculate the conditional Granger causality: -% -% det(restricted_residuals_cov_matrix) -% gc = ---------------------------------------------------- -% det(full_cov_matrix) -% -% see Barret, Barnett and Seth "Multivariate Granger causality and generalized variance". -% -% Call: -% connectivity = bst_granger_conditional(X, Y, 5, inputs); -% connectivity = bst_granger_conditional(X, [], 5, inputs); % every pair in X -% connectivity = bst_granger_conditional(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 - -% @============================================================================= -% This function is part of the Brainstorm software: -% https://neuroimage.usc.edu/brainstorm -% -% Copyright (c)2000-2020 University of Southern California & McGill University -% This software is distributed under the terms of the GNU General Public License -% as published by the Free Software Foundation. Further details on the GPLv3 -% license can be found at http://www.gnu.org/copyleft/gpl.html. -% -% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE -% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY -% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF -% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY -% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. -% -% For more information type "brainstorm license" at command prompt. -% =============================================================================@ -% -% Authors: Sergul Aydore & Syed Ashrafulla, 2012 -% Modified by: Davide Nuzzi, 2021 - -% default: 1 trial -if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) - inputs.nTrials = 1; -end - -% default: do not optimize order in MVAR modeling -if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) - inputs.flagFPE = false; -end - -% dimensions of the signals -nX = size(X, 1); -if ndims(X) == 3 - inputs.nTrials = size(X,3); - X = reshape(X, nX, []); -end -nSamples = size(X, 2); -nTimes = nSamples / inputs.nTrials; - -% if we are doing auto-causality, empty Y so that we only look at X -if ~isempty(Y) - nY = size(Y, 1); - Y = reshape(Y, nY, []); -end - -% remove linear effects if desired -if isfield(inputs, 'standardize') && inputs.standardize - detrender = [ ... - ones(1, nTimes); ... % constant trend in data - linspace(-1, 1, nTimes); ... % linear trend in data - 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data - ]; - - % detrend X - for iTrial = 1:inputs.nTrials - X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); - end - - % detrend Y only if it is not empty - if ~isempty(Y) - for iTrial = 1:inputs.nTrials - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); - end - end -end - -%% Iterate over all pairs of sinks & sources - -if isempty(Y) % auto-causality - - % setup - connectivity = zeros(nX); - - % the full model is evaluated only one time - - % multivariate model estimation - [transfers, noiseCovariance, order] = bst_mvar(X, order, inputs.nTrials, inputs.flagFPE); - - % data correlations using Yule-Walker (up to high order 50) - R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - - % only iterate over one triangle - for iX = 1:nX - - % iterate over all the pairs after iX - for iY = (iX+1):nX - - % restricted model iY -> iX - - % mask for the coefficients of the restricted model - mask = ones(size(X,1)); - mask(iX,iY) = 0; - - % restricted multivariate model (using masked row-by-row solution of YW equations) - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); - - % connectivity - connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); - - % restricted model iX -> iY - - % mask for the coefficients of the restricted model - mask = ones(size(X,1)); - mask(iY,iX) = 0; - - % restricted bivariate model (using masked row-by-row solution of - % YW equations) - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_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 - connectivity(iX, iX) = max([connectivity(iX, :) connectivity(:, iX)']); - - end - -else % cross-causality - - % setup - connectivity = zeros(nX, nY); - duplicates = zeros(0, 2); - - for iX = 1:nX - for iY = 1:nY - % X(iX,:) = sink, Y(iY,:) = source, Y(~iY,:) = conditional - % I have to check that the sink is different form all the other - % variables in Y - sink = X(iX,:); - source = Y(iY,:); - other_vars = Y(setdiff( 1:size(Y, 1), iY),:); - - % If sink and source are duplicates there's a correction - if max(abs(sink - source)) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere - - % If the target is equal to one of the conditioning - % variables, remove it - remove_inds = []; - for k = 1:size(other_vars,1) - if max(abs(sink - other_vars(k,:))) < eps - remove_inds = [remove_inds, k]; - end - end - other_vars(remove_inds,:) = []; - - - % multivariate model estimation (one sink, all sources) - [transfers, noiseCovariance, order] = bst_mvar([sink; source; other_vars], order, inputs.nTrials, inputs.flagFPE); - - % data correlations using Yule-Walker (up to high order 50) - R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - - % mask for the coefficients of the restricted model - mask = ones(2 + size(other_vars,1)); % size(other vars) + one source + one sink - mask(1,2) = 0; % Cut only the source -> sink coupling - - % restricted multivariate model (using masked row-by-row solution of YW equations) - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); - - % connectivity - connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); - - else % save duplicates to modify later - - duplicates(end+1, :) = [iX iY]; %#ok - - end - - end - end - - % for duplicate indices, set the causality value to the maximum over all inflows for iX and outflows for iY - for iDuplicate = 1:size(duplicates, 1) - connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2)) = ... - max([connectivity(duplicates(iDuplicate, 1), :) connectivity(:, duplicates(iDuplicate, 2))']); - end - -end - -%% Statistics: parametric p-values for causality in mean (based on regression coefficients) and variance (based on Wald statistics) - -% F statistic of connectivity when multiplied by number of regressors -pValue = 1; -% 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) - \ No newline at end of file From 727b0483130554afe75845e41ae4f45fbe1c02d5 Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 18:59:13 +0200 Subject: [PATCH 5/9] Delete bst_granger.m --- toolbox/connectivity/bst_granger.m | 218 ----------------------------- 1 file changed, 218 deletions(-) delete mode 100644 toolbox/connectivity/bst_granger.m diff --git a/toolbox/connectivity/bst_granger.m b/toolbox/connectivity/bst_granger.m deleted file mode 100644 index d8393798e..000000000 --- a/toolbox/connectivity/bst_granger.m +++ /dev/null @@ -1,218 +0,0 @@ -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 -% [X: MX x N or MX x N x T matrix] -% sources - second set of signals, one signal per row -% [Y: MY x N or MY x N x T matrix] -% (default: Y = X) -% order - maximum lag in AR model for causality in mean -% [p: nonnegative integer] -% inputs - structure of parameters: -% |-nTrials - # of trials in concantenated signal -% | [T: positive integer] -% |-standardize - if true (default), remove mean from each signal. -% | if false, assume signal has already been detrended -% |-flagFPE - if true, optimize order for AR model -% | if false (default), force same order in all AR models -% | [E: default false] -% -% Outputs: -% connectivity - A x B matrix of causalities from source to sink -% [C: MX x MY matrix] -% pValue - parametric p-value for corresponding Granger causality estimate -% [P: MX x MY matrix] -% -% 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: -% 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 - -% @============================================================================= -% This function is part of the Brainstorm software: -% https://neuroimage.usc.edu/brainstorm -% -% Copyright (c)2000-2020 University of Southern California & McGill University -% This software is distributed under the terms of the GNU General Public License -% as published by the Free Software Foundation. Further details on the GPLv3 -% license can be found at http://www.gnu.org/copyleft/gpl.html. -% -% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE -% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY -% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF -% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY -% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. -% -% For more information type "brainstorm license" at command prompt. -% =============================================================================@ -% -% Authors: Sergul Aydore & Syed Ashrafulla, 2012 -% Modified by: Davide Nuzzi, 2021 - -% default: 1 trial -if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) - inputs.nTrials = 1; -end - -% default: do not optimize order in MVAR modeling -if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) - inputs.flagFPE = false; -end - -% dimensions of the signals -nX = size(X, 1); -if ndims(X) == 3 - inputs.nTrials = size(X,3); - X = reshape(X, nX, []); -end -nSamples = size(X, 2); -nTimes = nSamples / inputs.nTrials; - -% if we are doing auto-causality, empty Y so that we only look at X -if ~isempty(Y) - nY = size(Y, 1); - Y = reshape(Y, nY, []); -end - -% remove linear effects if desired -if isfield(inputs, 'standardize') && inputs.standardize - detrender = [ ... - ones(1, nTimes); ... % constant trend in data - linspace(-1, 1, nTimes); ... % linear trend in data - 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data - ]; - - % detrend X - for iTrial = 1:inputs.nTrials - X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); - end - - % detrend Y only if it is not empty - if ~isempty(Y) - for iTrial = 1:inputs.nTrials - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); - end - end -end - -%% Iterate over all pairs of sinks & sources - -if isempty(Y) % auto-causality - - % setup - connectivity = zeros(nX); - - % only iterate over one triangle - for iX = 1:nX - - % iterate over all the pairs after iX - for iY = (iX+1):nX - - % 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 = YuleWalker_Inverse(transfers, noiseCovariance, 50); - - % restricted model iY -> iX - - % 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] = YuleWalker_Mask(R, mask); - - % connectivity - connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); - - % restricted model iX -> iY - - % 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] = YuleWalker_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 - connectivity(iX, iX) = max([connectivity(iX, :) connectivity(:, iX)']); - - end - -else % cross-causality - - % setup - connectivity = 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 - - % 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 = YuleWalker_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] = YuleWalker_Mask(R, mask); - - % connectivity - connectivity(iX, iY) = log(det(noiseCovariance_restricted) ./ det(noiseCovariance)); - - else % save duplicates to modify later - - duplicates(end+1, :) = [iX iY]; %#ok - - end - - end - end - - % for duplicate indices, set the causality value to the maximum over all inflows for iX and outflows for iY - for iDuplicate = 1:size(duplicates, 1) - connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2)) = ... - max([connectivity(duplicates(iDuplicate, 1), :) connectivity(:, duplicates(iDuplicate, 2))']); - end - -end - -%% Statistics: parametric p-values for causality in mean (based on regression coefficients) and variance (based on Wald statistics) - -% 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) - \ No newline at end of file From c98de7d3c1ac5cfa6e53f2a7555d008880eb3126 Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 18:59:21 +0200 Subject: [PATCH 6/9] Delete bst_granger_spectral.m --- toolbox/connectivity/bst_granger_spectral.m | 324 -------------------- 1 file changed, 324 deletions(-) delete mode 100644 toolbox/connectivity/bst_granger_spectral.m diff --git a/toolbox/connectivity/bst_granger_spectral.m b/toolbox/connectivity/bst_granger_spectral.m deleted file mode 100644 index b07ed5c3f..000000000 --- a/toolbox/connectivity/bst_granger_spectral.m +++ /dev/null @@ -1,324 +0,0 @@ -function [connectivity, pValues, freq] = bst_granger_spectral(X, Y, Fs, order, inputs) -% BST_GRANGER_SPECTRAL Granger causality at each frequency between any two -% signals. -% -% Inputs: -% X - first set of signals, one signal per row -% [X: A x N or A x N x T matrix] -% Y - second set of signals, one signal per row -% [Y: B x N or B x N x T matrix] -% (default: Y = X) -% Fs - sampling rate (we assume uniform sampling rate) -% [FS: scalar, FS > freq(end)*2] -% order - maximum order of autogressive model -% [p: integer > 1, default 10] -% inputs - structure of parameters: -% |-freq - frequencies of interest if desired -% |-freqResolution - maximum freq resolution in Hz, to limit NFFT -% | [DF: double, default [] (i.e. no limit)] -% |-nTrials - # of trials in concantenated signal -% |-flagFPE - if true, optimize order for autoregression -% | if false (default), use same order in autoregression -% |-standardize - if true (default), remove mean from each signal -% | if false, assume signal has already been detrended -% -% Outputs: -% connectivity - A x B matrix of spectral Granger causalities from -% source to sink. -% By default, GC(a,a) = 0 if Y is empty. -% [C: MX x MY x NF matrix] -% pValues - parametric p-value for corresponding spectral Granger -% causality in mean estimate (not implemented!) -% [P: MX x MY x NF matrix] -% freq - frequencies corresponding to the previous two metrics -% [F: NF x 1 vector] -% -% See also BST_GRANGER. -% -% Spectral causality measures are evaluated from the spectra of the full -% and restricted models as: -% -% gc(w) = det(S_restricted(w))/det(S_full(w)) -% -% see Cohen, Dror, et al. "A general spectral decomposition of causal influences -% applied to integrated information." Journal of neuroscience methods 330 (2020): 108443 -% for additional information. -% -% Call: -% connectivity = bst_granger_spectral(X, Y, 200, 10, inputs); % general call -% connectivity = bst_granger_spectral(X, [], 200, 30, inputs); % every pair -% Parameter examples: -% inputs.freq = 0:0.1:100; % specify desired frequencies -% inputs.freqResolution = 0.1; % have a high-point FFT -% inputs.nTrials = 9; % use trial-average covariances in AR estimation -% inputs.flagFPE = true; % use AR model with best information criteria - -% @============================================================================= -% This function is part of the Brainstorm software: -% https://neuroimage.usc.edu/brainstorm -% -% Copyright (c)2000-2020 University of Southern California & McGill University -% This software is distributed under the terms of the GNU General Public License -% as published by the Free Software Foundation. Further details on the GPLv3 -% license can be found at http://www.gnu.org/copyleft/gpl.html. -% -% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE -% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY -% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF -% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY -% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. -% -% For more information type "brainstorm license" at command prompt. -% =============================================================================@ -% -% Authors: Sergul Aydore & Syed Ashrafulla, 2012 -% Modified by: Davide Nuzzi, 2021 - -% reformat to a 2D matrix if necessary, and pull out # of trials -if ndims(X) == 3 - inputs.nTrials = size(X,3); - X = reshape(X, size(X, 1), []); - Y = reshape(Y, size(Y, 1), []); -elseif ~isfield(inputs, 'nTrials') - inputs.nTrials = 1; -end - -% lengths of things -nSamples = size(X, 2); -nTimes = nSamples / inputs.nTrials; - -% standardization: zero mean & unit variance, remove linear & quadratic trends as well -if isfield(inputs, 'standardize') && inputs.standardize - detrender = [ ... - ones(1, nTimes); ... % constant trend in data - linspace(-1, 1, nTimes); ... % linear trend in data - 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data - ]; - - % detrend X - for iTrial = 1:inputs.nTrials - X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); - end - - % detrend Y only if it is not empty - if ~isempty(Y) - for iTrial = 1:inputs.nTrials - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); - end - end -end - -% number of FFT bins required -if ~isfield(inputs, 'freq') || isempty(inputs.freq) % frequencies of interest are not defined - if isfield(inputs, 'freqResolution') && ~isempty(inputs.freqResolution) && (size(X,2) > round(Fs / inputs.freqResolution)) - nFFT = 2^nextpow2( round(Fs / inputs.freqResolution) ); - else % use a default frequency resolution of 1Hz to mirror the standard frequency resolution in bst_coherence_welch.m - nFFT = 2^nextpow2( round(Fs / 1) ); - end -elseif numel(inputs.freq) == 1 - nFFT = inputs.freq; -else - nFFT = 2^nextpow2(max( length(inputs.freq)-1, (Fs/2) / min(diff( sort(inputs.freq(:)) )) )) * 2; -end - -% default: Order 10 used in BrainStorm -if isempty(order) - order = 10; -end - -% default: Single-trial -if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) - inputs.nTrials = 1; -end - -% default: do not optimize model order -if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) - inputs.flagFPE = false; -end - -%% Differentiate between auto-causality and cross-causality for speed - -if isempty(Y) % auto-causality between signals in X, so we can halve the number of models to estimate - - % preallocate the spectral causality matrix - connectivity = zeros(size(X, 1), size(X, 1), nFFT/2); - - % iterate over all pairs of sinks & sources - for iX = 1:size(X, 1) - for iY = iX+1 : size(X, 1) % to avoid auto-causality - - % two-variate model - [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); X(iY, :)], order, inputs.nTrials, inputs.flagFPE); - - % spectra of full model - [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); - - % data correlations using Yule-Walker (up to high order 50) - R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - - % restricted model iY -> iX - - % 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] = YuleWalker_Mask(R, mask); - - % spectra of restricted system - [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); - - % connectivity at each frequence (see Cohen et al., 2020) - for n = 1:length(freq) - connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger - end - - % restricted model iX -> iY - - % mask for the coefficients of the restricted model - mask = ones(2); - mask(2,1) = 0; - - % restricted bivariate model - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); - - % spectra of restricted system - [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); - - % connectivity at each frequence (see Cohen et al., 2020) - for n = 1:length(freq) - connectivity(iY, iX, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger - end - end - - % diagonal will equal the maximum of all inflows and outflows for iX, specific to each frequency - connectivity(iX, iX, :) = max( max(connectivity(iX, :, :), [], 2), max(connectivity(:, iX, :), [], 1) ); - - end - -else % we have to use all pairs of signals - - % preallocate the spectral causality matrix - connectivity = zeros(size(X, 1), size(Y, 1), nFFT/2); - duplicates = zeros(0, 2); - - % iterate over all pairs of sinks & sources - for iX = 1:size(X, 1) - for iY = 1:size(Y, 1) - - if max(abs(X(iX, :) - Y(iY, :))) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere - - % two-variate model - [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); Y(iY, :)], order, inputs.nTrials, inputs.flagFPE); - - % spectra of full model - [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); - - % data correlations using Yule-Walker (up to high order 50) - R = YuleWalker_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] = YuleWalker_Mask(R, mask); - - % spectra of restricted system - [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); - - % connectivity at each frequence (see Cohen et al., 2020) - for n = 1:length(freq) - connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger - end - - else % save duplicates to modify later - - duplicates(end+1, :) = [iX iY]; %#ok - - end - - end - end - - % for duplicate indices, set the causality value to the maximum of all inflows for iX and outflows for iY - for iDuplicate = 1:size(duplicates, 1) - connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2), :) = max( ... - max(connectivity(duplicates(iDuplicate, 1), :, :), [], 2), ... - max(connectivity(:, duplicates(iDuplicate, 2), :), [], 1) ... - ); - end - - -end - -%% Interpolate to desired frequencies & perform statistics if desired -if isfield(inputs, 'freq') && ~isempty(inputs.freq) && numel(inputs.freq) > 1 - connectivity = permute(interp1(freq, permute(connectivity, [3 1 2]), inputs.freq), [2 3 1]); - % pValues = permute(interp1(freq, permute(pValues, [3 1 2]), inputs.freq), [2 3 1]); - freq = inputs.freq; -end - -pValues = NaN; % no parametric p-values for now - -end - -%% ======================================================== spectra estimation ======================================================== -function [spectra, freq] = bst_granger_spectral_spectrum(A, Sigma, nFFT, Fs) -% BST_GRANGER_SPECTRAL_SPECTRUM Calculate the parametric spectra of a multivariate system -% with given MVAR coefficients and an estimate of the -% covariance matrix of the innovation (i.e. noise) -% process. -% -% Inputs: -% transfers - transfer matrices in AR process -% [A: N x NP matrix, P = order] -% noiseCovariance - variance of residuals -% [C: N x N matrix] -% nFFT - number of FFT bins to calculate spectra -% [NF: positive number, usually power of 2] -% Fs - sampling frequency of data -% [FS: double, FS > freq(end)*2] -% -% Outputs: -% spectra - cross-spectrum between each pair of variables -% [S: N x N x NF/2 matrix] -% freq - frequencies used based on # of FFT bins -% [F: NF/2 x 1 matrix] -% --> all outputs have length NF/2, ignoring frequenices past Fs/2 -% Call: -% spectra = bst_mvar_spectrum(transfers, C, 512, 200); % basic -% spectra = bst_mvar_spectrum(transfers, C, 2048, 200); % add FFT interp - -% frequencies to estimate cross-spectra -freq = Fs/2*linspace(0, 1, nFFT/2 + 1); -freq(end) = []; - -% get number of variables, fft bins and order of the model -M = nFFT/2; -[N,pN] = size(A); -p = pN/N; - -% evaluate the trasfer function and the spectra -H = complex(zeros(N,N,M)); -spectra = complex(zeros(N,N,M)); - -A_r = reshape(A,[N,N,p]); -w_vec = 0:pi/(nFFT/2-1):pi; - -for n = 1:M - e = zeros(p,1); - for k = 1:p - e(k) = exp(-1i * w_vec(n) * k); - end - e = permute(repmat(e,[1,N,N]),[2,3,1]); - A_w = eye(N) - sum(A_r .* e,3); - H(:,:,n) = inv(A_w); - spectra(:,:,n) = H(:,:,n) * Sigma * ctranspose(H(:,:,n)); -end - - -end %% <== FUNCTION END \ No newline at end of file From 7373c86819f407c55ee8d4505dc017f6e4d8d9fe Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 18:59:28 +0200 Subject: [PATCH 7/9] Delete bst_granger_spectral_conditional.m --- .../bst_granger_spectral_conditional.m | 350 ------------------ 1 file changed, 350 deletions(-) delete mode 100644 toolbox/connectivity/bst_granger_spectral_conditional.m diff --git a/toolbox/connectivity/bst_granger_spectral_conditional.m b/toolbox/connectivity/bst_granger_spectral_conditional.m deleted file mode 100644 index f2015044b..000000000 --- a/toolbox/connectivity/bst_granger_spectral_conditional.m +++ /dev/null @@ -1,350 +0,0 @@ -function [connectivity, pValues, freq] = bst_granger_spectral_conditional(X, Y, Fs, order, inputs) -% BST_GRANGER_SPECTRAL_CONDITIONAL Granger causality at each frequency between any two -% signals, conditioned on all the remaining signals -% -% Inputs: -% X - first set of signals, one signal per row -% [X: A x N or A x N x T matrix] -% Y - second set of signals, one signal per row -% [Y: B x N or B x N x T matrix] -% (default: Y = X) -% Fs - sampling rate (we assume uniform sampling rate) -% [FS: scalar, FS > freq(end)*2] -% order - maximum order of autogressive model -% [p: integer > 1, default 10] -% inputs - structure of parameters: -% |-freq - frequencies of interest if desired -% |-freqResolution - maximum freq resolution in Hz, to limit NFFT -% | [DF: double, default [] (i.e. no limit)] -% |-nTrials - # of trials in concantenated signal -% |-flagFPE - if true, optimize order for autoregression -% | if false (default), use same order in autoregression -% |-standardize - if true (default), remove mean from each signal -% | if false, assume signal has already been detrended -% -% Outputs: -% connectivity - A x B matrix of spectral Granger causalities from -% source to sink, conditioned on all the other variables. -% If Y is empty, the matrix contains the spectral GC -% from each source variable to each sink variable, -% conditioned on the remaining variable. -% If Y is not empty the matrix contains the spectral -% causalities from each variable in Y to each -% variable in X, conditionend on the other variables -% in Y. -% [C: MX x MY x NF matrix] -% pValues - parametric p-value for corresponding spectral Granger -% causality in mean estimate (not implemented yet!) -% [P: MX x MY x NF matrix] -% freq - frequencies corresponding to the previous two metrics -% [F: NF x 1 vector] -% -% See also BST_GRANGER_CONDITIONAL. -% - -% Spectral causality measures are evaluated from the spectra of the full -% and restricted models as: -% -% gc(w) = det(S_restricted(w))/det(S_full(w)) -% -% see Cohen, Dror, et al. "A general spectral decomposition of causal influences -% applied to integrated information." Journal of neuroscience methods 330 (2020): 108443 -% for additional information. -% -% Call: -% connectivity = bst_granger_spectral_conditional(X, Y, 200, 10, inputs); % general call -% connectivity = bst_granger_spectral_conditional(X, [], 200, 30, inputs); % every pair in X -% Parameter examples: -% inputs.freq = 0:0.1:100; % specify desired frequencies -% inputs.freqResolution = 0.1; % have a high-point FFT -% inputs.nTrials = 9; % use trial-average covariances in AR estimation -% inputs.flagFPE = true; % use AR model with best information criteria - -% @============================================================================= -% This function is part of the Brainstorm software: -% https://neuroimage.usc.edu/brainstorm -% -% Copyright (c)2000-2020 University of Southern California & McGill University -% This software is distributed under the terms of the GNU General Public License -% as published by the Free Software Foundation. Further details on the GPLv3 -% license can be found at http://www.gnu.org/copyleft/gpl.html. -% -% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE -% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY -% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF -% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY -% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. -% -% For more information type "brainstorm license" at command prompt. -% =============================================================================@ -% -% Authors: Sergul Aydore & Syed Ashrafulla, 2012 -% Modified by: Davide Nuzzi, 2021 - -% reformat to a 2D matrix if necessary, and pull out # of trials -if ndims(X) == 3 - inputs.nTrials = size(X,3); - X = reshape(X, size(X, 1), []); - Y = reshape(Y, size(Y, 1), []); -elseif ~isfield(inputs, 'nTrials') - inputs.nTrials = 1; -end - -% lengths of things -nSamples = size(X, 2); -nTimes = nSamples / inputs.nTrials; - -% standardization: zero mean & unit variance, remove linear & quadratic trends as well -if isfield(inputs, 'standardize') && inputs.standardize - detrender = [ ... - ones(1, nTimes); ... % constant trend in data - linspace(-1, 1, nTimes); ... % linear trend in data - 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data - ]; - - % detrend X - for iTrial = 1:inputs.nTrials - X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); - end - - % detrend Y only if it is not empty - if ~isempty(Y) - for iTrial = 1:inputs.nTrials - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; - Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); - end - end -end - -% number of FFT bins required -if ~isfield(inputs, 'freq') || isempty(inputs.freq) % frequencies of interest are not defined - if isfield(inputs, 'freqResolution') && ~isempty(inputs.freqResolution) && (size(X,2) > round(Fs / inputs.freqResolution)) - nFFT = 2^nextpow2( round(Fs / inputs.freqResolution) ); - else % use a default frequency resolution of 1Hz to mirror the standard frequency resolution in bst_coherence_welch.m - nFFT = 2^nextpow2( round(Fs / 1) ); - end -elseif numel(inputs.freq) == 1 - nFFT = inputs.freq; -else - nFFT = 2^nextpow2(max( length(inputs.freq)-1, (Fs/2) / min(diff( sort(inputs.freq(:)) )) )) * 2; -end - -% default: Order 10 used in BrainStorm -if isempty(order) - order = 10; -end - -% default: Single-trial -if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) - inputs.nTrials = 1; -end - -% default: do not optimize model order -if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) - inputs.flagFPE = false; -end - -%% Differentiate between auto-causality and cross-causality for speed - -if isempty(Y) % causality between signals in X - - % preallocate the spectral causality matrix - connectivity = zeros(size(X, 1), size(X, 1), nFFT/2); - - % the full model is evaluated only one time - - % multivariate model estimation - [transfers, noiseCovariance, order] = bst_mvar(X, order, inputs.nTrials, inputs.flagFPE); - - % spectra for the full model - [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); - - % data correlations using Yule-Walker (up to high order 50) - R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - - % iterate over all pairs of sinks & sources - for iX = 1:size(X, 1) - for iY = iX+1 : size(X, 1) % to avoid auto-causality - - % restricted model iY -> iX - - % mask for the coefficients of the restricted model - mask = ones(size(X,1)); - mask(iX,iY) = 0; - - % restricted multivariate model (using masked row-by-row solution of YW equations) - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); - - % spectra of restricted system - [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); - - % connectivity at each frequence (see Cohen et al., 2020) - for n = 1:length(freq) - connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger - end - - % restricted model iX -> iY - - % mask for the coefficients of the restricted model - mask = ones(size(X,1)); - mask(iY,iX) = 0; - - % restricted multivariate model (using masked row-by-row solution of YW equations) - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); - - % spectra of restricted system - [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); - - % connectivity at each frequence (see Cohen et al., 2020) - for n = 1:length(freq) - connectivity(iY, iX, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger - end - end - - % diagonal will equal the maximum of all inflows and outflows for iX, specific to each frequency - connectivity(iX, iX, :) = max( max(connectivity(iX, :, :), [], 2), max(connectivity(:, iX, :), [], 1) ); - end - -else % we have to use all pairs of signals - - % preallocate the spectral causality matrix - connectivity = zeros(size(X, 1), size(Y, 1), nFFT/2); - duplicates = zeros(0, 2); - - % iterate over all pairs of sinks & sources - for iX = 1:size(X, 1) - for iY = 1:size(Y, 1) - - % X(iX,:) = sink, Y(iY,:) = source, Y(~iY,:) = conditional - % I have to check that the sink is different form all the other - % variables in Y - sink = X(iX,:); - source = Y(iY,:); - other_vars = Y(setdiff( 1:size(Y, 1), iY),:); - - % If sink and source are duplicates there's a correction - % - if max(abs(sink - source)) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere - - % If the target is equal to one of the conditioning - % variables, remove it - remove_inds = []; - for k = 1:size(other_vars,1) - if max(abs(sink - other_vars(k,:))) < eps - remove_inds = [remove_inds, k]; - end - end - other_vars(remove_inds,:) = []; - - - % multivariate model estimation (one sink, all sources) - [transfers, noiseCovariance, order] = bst_mvar([sink; source; other_vars], order, inputs.nTrials, inputs.flagFPE); - - % spectra for the full model - [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); - - % data correlations using Yule-Walker (up to high order 50) - R = YuleWalker_Inverse(transfers, noiseCovariance, 50); - - % mask for the coefficients of the restricted model - mask = ones(2 + size(other_vars,1)); % size(other vars) + one source + one sink - mask(1,2) = 0; % Cut only the source -> sink coupling - - % restricted multivariate model (using masked row-by-row solution of YW equations) - [transfers_restricted,noiseCovariance_restricted] = YuleWalker_Mask(R, mask); - - % spectra of restricted system - [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); - - % connectivity at each frequence (see Cohen et al., 2020) - for n = 1:length(freq) - connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger - end - - else % save duplicates to modify later - - duplicates(end+1, :) = [iX iY]; %#ok - - end - - end - end - - % for duplicate indices, set the causality value to the maximum of all inflows for iX and outflows for iY - for iDuplicate = 1:size(duplicates, 1) - connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2), :) = max( ... - max(connectivity(duplicates(iDuplicate, 1), :, :), [], 2), ... - max(connectivity(:, duplicates(iDuplicate, 2), :), [], 1) ... - ); - end - -end - -%% Interpolate to desired frequencies & perform statistics if desired -if isfield(inputs, 'freq') && ~isempty(inputs.freq) && numel(inputs.freq) > 1 - connectivity = permute(interp1(freq, permute(connectivity, [3 1 2]), inputs.freq), [2 3 1]); - % pValues = permute(interp1(freq, permute(pValues, [3 1 2]), inputs.freq), [2 3 1]); - freq = inputs.freq; -end - -pValues = NaN; % no parametric p-values for now - -end - -%% ======================================================== spectra estimation ======================================================== -function [spectra, freq] = bst_granger_spectral_spectrum(A, Sigma, nFFT, Fs) -% BST_GRANGER_SPECTRAL_SPECTRUM Calculate the parametric spectra of a multivariate system -% with given MVAR coefficients and an estimate of the -% covariance matrix of the innovation (i.e. noise) -% process. -% -% Inputs: -% transfers - transfer matrices in AR process -% [A: N x NP matrix, P = order] -% noiseCovariance - variance of residuals -% [C: N x N matrix] -% nFFT - number of FFT bins to calculate spectra -% [NF: positive number, usually power of 2] -% Fs - sampling frequency of data -% [FS: double, FS > freq(end)*2] -% -% Outputs: -% spectra - cross-spectrum between each pair of variables -% [S: N x N x NF/2 matrix] -% freq - frequencies used based on # of FFT bins -% [F: NF/2 x 1 matrix] -% --> all outputs have length NF/2, ignoring frequenices past Fs/2 -% Call: -% spectra = bst_mvar_spectrum(transfers, C, 512, 200); % basic -% spectra = bst_mvar_spectrum(transfers, C, 2048, 200); % add FFT interp - -% frequencies to estimate cross-spectra -freq = Fs/2*linspace(0, 1, nFFT/2 + 1); -freq(end) = []; - -% get number of variables, fft bins and order of the model -M = nFFT/2; -[N,pN] = size(A); -p = pN/N; - -% evaluate the trasfer function and the spectra -H = complex(zeros(N,N,M)); -spectra = complex(zeros(N,N,M)); - -A_r = reshape(A,[N,N,p]); -w_vec = 0:pi/(nFFT/2-1):pi; - -for n = 1:M - e = zeros(p,1); - for k = 1:p - e(k) = exp(-1i * w_vec(n) * k); - end - e = permute(repmat(e,[1,N,N]),[2,3,1]); - A_w = eye(N) - sum(A_r .* e,3); - H(:,:,n) = inv(A_w); - spectra(:,:,n) = H(:,:,n) * Sigma * ctranspose(H(:,:,n)); -end - - -end %% <== FUNCTION END \ No newline at end of file From e0014f37c9772cce4b095c5ca220e76e8e9c42ea Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 19:00:12 +0200 Subject: [PATCH 8/9] Add files via upload --- toolbox/connectivity/bst_granger.m | 218 +++++++++++++ toolbox/connectivity/bst_granger_spectral.m | 324 ++++++++++++++++++++ 2 files changed, 542 insertions(+) create mode 100644 toolbox/connectivity/bst_granger.m create mode 100644 toolbox/connectivity/bst_granger_spectral.m diff --git a/toolbox/connectivity/bst_granger.m b/toolbox/connectivity/bst_granger.m new file mode 100644 index 000000000..653d4fd1c --- /dev/null +++ b/toolbox/connectivity/bst_granger.m @@ -0,0 +1,218 @@ +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 +% [X: MX x N or MX x N x T matrix] +% sources - second set of signals, one signal per row +% [Y: MY x N or MY x N x T matrix] +% (default: Y = X) +% order - maximum lag in AR model for causality in mean +% [p: nonnegative integer] +% inputs - structure of parameters: +% |-nTrials - # of trials in concantenated signal +% | [T: positive integer] +% |-standardize - if true (default), remove mean from each signal. +% | if false, assume signal has already been detrended +% |-flagFPE - if true, optimize order for AR model +% | if false (default), force same order in all AR models +% | [E: default false] +% +% Outputs: +% connectivity - A x B matrix of causalities from source to sink +% [C: MX x MY matrix] +% pValue - parametric p-value for corresponding Granger causality estimate +% [P: MX x MY matrix] +% +% 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: +% 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 + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Sergul Aydore & Syed Ashrafulla, 2012 +% Modified by: Davide Nuzzi, 2021 + +% default: 1 trial +if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) + inputs.nTrials = 1; +end + +% default: do not optimize order in MVAR modeling +if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) + inputs.flagFPE = false; +end + +% dimensions of the signals +nX = size(X, 1); +if ndims(X) == 3 + inputs.nTrials = size(X,3); + X = reshape(X, nX, []); +end +nSamples = size(X, 2); +nTimes = nSamples / inputs.nTrials; + +% if we are doing auto-causality, empty Y so that we only look at X +if ~isempty(Y) + nY = size(Y, 1); + Y = reshape(Y, nY, []); +end + +% remove linear effects if desired +if isfield(inputs, 'standardize') && inputs.standardize + detrender = [ ... + ones(1, nTimes); ... % constant trend in data + linspace(-1, 1, nTimes); ... % linear trend in data + 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data + ]; + + % detrend X + for iTrial = 1:inputs.nTrials + X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); + end + + % detrend Y only if it is not empty + if ~isempty(Y) + for iTrial = 1:inputs.nTrials + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); + end + end +end + +%% Iterate over all pairs of sinks & sources + +if isempty(Y) % auto-causality + + % setup + connectivity = zeros(nX); + + % only iterate over one triangle + for iX = 1:nX + + % iterate over all the pairs after iX + for iY = (iX+1):nX + + % 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 + + % 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)); + + % restricted model iX -> iY + + % 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 + connectivity(iX, iX) = max([connectivity(iX, :) connectivity(:, iX)']); + + end + +else % cross-causality + + % setup + connectivity = 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 + + % 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 + + duplicates(end+1, :) = [iX iY]; %#ok + + end + + end + end + + % for duplicate indices, set the causality value to the maximum over all inflows for iX and outflows for iY + for iDuplicate = 1:size(duplicates, 1) + connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2)) = ... + max([connectivity(duplicates(iDuplicate, 1), :) connectivity(:, duplicates(iDuplicate, 2))']); + end + +end + +%% Statistics: parametric p-values for causality in mean (based on regression coefficients) and variance (based on Wald statistics) + +% 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) + \ No newline at end of file diff --git a/toolbox/connectivity/bst_granger_spectral.m b/toolbox/connectivity/bst_granger_spectral.m new file mode 100644 index 000000000..f3605eb7f --- /dev/null +++ b/toolbox/connectivity/bst_granger_spectral.m @@ -0,0 +1,324 @@ +function [connectivity, pValues, freq] = bst_granger_spectral(X, Y, Fs, order, inputs) +% BST_GRANGER_SPECTRAL Granger causality at each frequency between any two +% signals. +% +% Inputs: +% X - first set of signals, one signal per row +% [X: A x N or A x N x T matrix] +% Y - second set of signals, one signal per row +% [Y: B x N or B x N x T matrix] +% (default: Y = X) +% Fs - sampling rate (we assume uniform sampling rate) +% [FS: scalar, FS > freq(end)*2] +% order - maximum order of autogressive model +% [p: integer > 1, default 10] +% inputs - structure of parameters: +% |-freq - frequencies of interest if desired +% |-freqResolution - maximum freq resolution in Hz, to limit NFFT +% | [DF: double, default [] (i.e. no limit)] +% |-nTrials - # of trials in concantenated signal +% |-flagFPE - if true, optimize order for autoregression +% | if false (default), use same order in autoregression +% |-standardize - if true (default), remove mean from each signal +% | if false, assume signal has already been detrended +% +% Outputs: +% connectivity - A x B matrix of spectral Granger causalities from +% source to sink. +% By default, GC(a,a) = 0 if Y is empty. +% [C: MX x MY x NF matrix] +% pValues - parametric p-value for corresponding spectral Granger +% causality in mean estimate (not implemented!) +% [P: MX x MY x NF matrix] +% freq - frequencies corresponding to the previous two metrics +% [F: NF x 1 vector] +% +% See also BST_GRANGER. +% +% Spectral causality measures are evaluated from the spectra of the full +% and restricted models as: +% +% gc(w) = det(S_restricted(w))/det(S_full(w)) +% +% see Cohen, Dror, et al. "A general spectral decomposition of causal influences +% applied to integrated information." Journal of neuroscience methods 330 (2020): 108443 +% for additional information. +% +% Call: +% connectivity = bst_granger_spectral(X, Y, 200, 10, inputs); % general call +% connectivity = bst_granger_spectral(X, [], 200, 30, inputs); % every pair +% Parameter examples: +% inputs.freq = 0:0.1:100; % specify desired frequencies +% inputs.freqResolution = 0.1; % have a high-point FFT +% inputs.nTrials = 9; % use trial-average covariances in AR estimation +% inputs.flagFPE = true; % use AR model with best information criteria + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Sergul Aydore & Syed Ashrafulla, 2012 +% Modified by: Davide Nuzzi, 2021 + +% reformat to a 2D matrix if necessary, and pull out # of trials +if ndims(X) == 3 + inputs.nTrials = size(X,3); + X = reshape(X, size(X, 1), []); + Y = reshape(Y, size(Y, 1), []); +elseif ~isfield(inputs, 'nTrials') + inputs.nTrials = 1; +end + +% lengths of things +nSamples = size(X, 2); +nTimes = nSamples / inputs.nTrials; + +% standardization: zero mean & unit variance, remove linear & quadratic trends as well +if isfield(inputs, 'standardize') && inputs.standardize + detrender = [ ... + ones(1, nTimes); ... % constant trend in data + linspace(-1, 1, nTimes); ... % linear trend in data + 3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data + ]; + + % detrend X + for iTrial = 1:inputs.nTrials + X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes)); + end + + % detrend Y only if it is not empty + if ~isempty(Y) + for iTrial = 1:inputs.nTrials + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender; + Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes)); + end + end +end + +% number of FFT bins required +if ~isfield(inputs, 'freq') || isempty(inputs.freq) % frequencies of interest are not defined + if isfield(inputs, 'freqResolution') && ~isempty(inputs.freqResolution) && (size(X,2) > round(Fs / inputs.freqResolution)) + nFFT = 2^nextpow2( round(Fs / inputs.freqResolution) ); + else % use a default frequency resolution of 1Hz to mirror the standard frequency resolution in bst_coherence_welch.m + nFFT = 2^nextpow2( round(Fs / 1) ); + end +elseif numel(inputs.freq) == 1 + nFFT = inputs.freq; +else + nFFT = 2^nextpow2(max( length(inputs.freq)-1, (Fs/2) / min(diff( sort(inputs.freq(:)) )) )) * 2; +end + +% default: Order 10 used in BrainStorm +if isempty(order) + order = 10; +end + +% default: Single-trial +if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials) + inputs.nTrials = 1; +end + +% default: do not optimize model order +if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE) + inputs.flagFPE = false; +end + +%% Differentiate between auto-causality and cross-causality for speed + +if isempty(Y) % auto-causality between signals in X, so we can halve the number of models to estimate + + % preallocate the spectral causality matrix + connectivity = zeros(size(X, 1), size(X, 1), nFFT/2); + + % iterate over all pairs of sinks & sources + for iX = 1:size(X, 1) + for iY = iX+1 : size(X, 1) % to avoid auto-causality + + % two-variate model + [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); X(iY, :)], order, inputs.nTrials, inputs.flagFPE); + + % spectra of full model + [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); + + % data correlations using Yule-Walker (up to high order 50) + R = yule_walker_inverse(transfers, noiseCovariance, 50); + + % restricted model iY -> iX + + % 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); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + + % restricted model iX -> iY + + % mask for the coefficients of the restricted model + mask = ones(2); + mask(2,1) = 0; + + % restricted bivariate model + [transfers_restricted,noiseCovariance_restricted] = yule_walker_mask(R, mask); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iY, iX, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + end + + % diagonal will equal the maximum of all inflows and outflows for iX, specific to each frequency + connectivity(iX, iX, :) = max( max(connectivity(iX, :, :), [], 2), max(connectivity(:, iX, :), [], 1) ); + + end + +else % we have to use all pairs of signals + + % preallocate the spectral causality matrix + connectivity = zeros(size(X, 1), size(Y, 1), nFFT/2); + duplicates = zeros(0, 2); + + % iterate over all pairs of sinks & sources + for iX = 1:size(X, 1) + for iY = 1:size(Y, 1) + + if max(abs(X(iX, :) - Y(iY, :))) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere + + % two-variate model + [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); Y(iY, :)], order, inputs.nTrials, inputs.flagFPE); + + % spectra of full model + [spectra, freq] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs); + + % 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); + + % spectra of restricted system + [spectra_restricted, ~] = bst_granger_spectral_spectrum(transfers_restricted, noiseCovariance_restricted, nFFT, Fs); + + % connectivity at each frequence (see Cohen et al., 2020) + for n = 1:length(freq) + connectivity(iX, iY, n) = log(abs(det(spectra_restricted(:,:,n))) ./ abs(det(spectra(:,:,n))));% Geweke-Granger + end + + else % save duplicates to modify later + + duplicates(end+1, :) = [iX iY]; %#ok + + end + + end + end + + % for duplicate indices, set the causality value to the maximum of all inflows for iX and outflows for iY + for iDuplicate = 1:size(duplicates, 1) + connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2), :) = max( ... + max(connectivity(duplicates(iDuplicate, 1), :, :), [], 2), ... + max(connectivity(:, duplicates(iDuplicate, 2), :), [], 1) ... + ); + end + + +end + +%% Interpolate to desired frequencies & perform statistics if desired +if isfield(inputs, 'freq') && ~isempty(inputs.freq) && numel(inputs.freq) > 1 + connectivity = permute(interp1(freq, permute(connectivity, [3 1 2]), inputs.freq), [2 3 1]); + % pValues = permute(interp1(freq, permute(pValues, [3 1 2]), inputs.freq), [2 3 1]); + freq = inputs.freq; +end + +pValues = NaN; % no parametric p-values for now + +end + +%% ======================================================== spectra estimation ======================================================== +function [spectra, freq] = bst_granger_spectral_spectrum(A, Sigma, nFFT, Fs) +% BST_GRANGER_SPECTRAL_SPECTRUM Calculate the parametric spectra of a multivariate system +% with given MVAR coefficients and an estimate of the +% covariance matrix of the innovation (i.e. noise) +% process. +% +% Inputs: +% transfers - transfer matrices in AR process +% [A: N x NP matrix, P = order] +% noiseCovariance - variance of residuals +% [C: N x N matrix] +% nFFT - number of FFT bins to calculate spectra +% [NF: positive number, usually power of 2] +% Fs - sampling frequency of data +% [FS: double, FS > freq(end)*2] +% +% Outputs: +% spectra - cross-spectrum between each pair of variables +% [S: N x N x NF/2 matrix] +% freq - frequencies used based on # of FFT bins +% [F: NF/2 x 1 matrix] +% --> all outputs have length NF/2, ignoring frequenices past Fs/2 +% Call: +% spectra = bst_mvar_spectrum(transfers, C, 512, 200); % basic +% spectra = bst_mvar_spectrum(transfers, C, 2048, 200); % add FFT interp + +% frequencies to estimate cross-spectra +freq = Fs/2*linspace(0, 1, nFFT/2 + 1); +freq(end) = []; + +% get number of variables, fft bins and order of the model +M = nFFT/2; +[N,pN] = size(A); +p = pN/N; + +% evaluate the trasfer function and the spectra +H = complex(zeros(N,N,M)); +spectra = complex(zeros(N,N,M)); + +A_r = reshape(A,[N,N,p]); +w_vec = 0:pi/(nFFT/2-1):pi; + +for n = 1:M + e = zeros(p,1); + for k = 1:p + e(k) = exp(-1i * w_vec(n) * k); + end + e = permute(repmat(e,[1,N,N]),[2,3,1]); + A_w = eye(N) - sum(A_r .* e,3); + H(:,:,n) = inv(A_w); + spectra(:,:,n) = H(:,:,n) * Sigma * ctranspose(H(:,:,n)); +end + + +end %% <== FUNCTION END \ No newline at end of file From 1e7ca9eae964f979e08fd96947ccc1b665ee7fbe Mon Sep 17 00:00:00 2001 From: DavideNuzzi <37365594+DavideNuzzi@users.noreply.github.com> Date: Wed, 18 Aug 2021 19:01:03 +0200 Subject: [PATCH 9/9] Add files via upload --- .../private/yule_walker_inverse.m | 43 +++++++++++++++ .../connectivity/private/yule_walker_mask.m | 52 +++++++++++++++++++ 2 files changed, 95 insertions(+) create mode 100644 toolbox/connectivity/private/yule_walker_inverse.m create mode 100644 toolbox/connectivity/private/yule_walker_mask.m diff --git a/toolbox/connectivity/private/yule_walker_inverse.m b/toolbox/connectivity/private/yule_walker_inverse.m new file mode 100644 index 000000000..575ab657a --- /dev/null +++ b/toolbox/connectivity/private/yule_walker_inverse.m @@ -0,0 +1,43 @@ +function R = yule_walker_inverse(A,Sigma, maxOrder) +% YULE_WALKER_INVErSE Solve the inverse Yule-Walker system of equations. +% Given the couplings A and the covariance matrix +% Sigma it evaluated the correlations of the time +% series up to an arbitrary order +% +% Inputs: +% A - Couplings of the VAR model +% [N x Np, where Np = N x p] +% Sigma - Residuals covariance matrix +% [N x N] +% maxOrder - Maximum order for the data correlations +% +% Outputs: +% R - Correlations of the time series obtained using the +% Yule-Walker equations +% [N x N x (maxOrder+1)] + + [N,~] = size(Sigma); + p = size(A,2)/N; + + Im = eye(N*p); + F = [A; Im(1:end-N,:)]; + Delta=zeros(N*p); + Delta(1:N,1:N) = Sigma; + BigSigma = dlyap(F,Delta); + + R=NaN*ones(N,N,maxOrder+1); + for i=1:p + R(:,:,i)=BigSigma(1:N,N*(i-1)+1:N*i); + end + + for k=p+1:maxOrder+1 + Rk=R(:,:,k-1:-1:k-p); + Rm=[]; + for ki=1:p + Rm=[Rm; Rk(:,:,ki)]; + end + R(:,:,k)=A*Rm; + end + +end + diff --git a/toolbox/connectivity/private/yule_walker_mask.m b/toolbox/connectivity/private/yule_walker_mask.m new file mode 100644 index 000000000..90be98a5a --- /dev/null +++ b/toolbox/connectivity/private/yule_walker_mask.m @@ -0,0 +1,52 @@ +function [A,Sigma] = yule_walker_mask(Gamma, A_mask) + +% YULE_WALKER_MASK Solve the Yule-Walker system of equations under +% some constraints on the couplings. +% +% Inputs: +% Gamma - Correlations of the time series +% [N x N x (p+1)] +% A_mask - Mask to cut off some couplings +% [N x N] +% +% Outputs: +% A - Couplings of the VAR model +% [N x Np, where Np = N x p] +% Sigma - Residuals covariance matrix +% [N x N] + +[N,~,q] = size(Gamma); +p = q - 1; + +% Initialize the matrices used for the evaluation +Psi = zeros(N*p); +G = reshape(Gamma(:,:,2:end),[N,N*p]); +A_mask = repmat(A_mask,[1,p]); + +for i = 1:p + for j = 1:p + k = j - i; + if (k > 0) % Colonna > Riga = metà superiore destra + Psi((i-1)*N+1:i*N, (j-1)*N+1:j*N) = Gamma(:,:,k + 1); + else + Psi((i-1)*N+1:i*N, (j-1)*N+1:j*N) = Gamma(:,:,abs(k) + 1)'; + end + end +end + +% Solve row-by-row +A = zeros(N,N*p); + +for i = 1:N + % Solve only for the indices that aren't masked off + ind = find(A_mask(i,:) ~= 0); + Psi_t = Psi(ind,ind); + G_t = G(i,ind); + + A(i,ind) = G_t / Psi_t; +end + +Sigma = Gamma(:,:,1) - A * G'; + +end +