diff --git a/toolbox/connectivity/bst_granger.m b/toolbox/connectivity/bst_granger.m index e3fa4b4d4..653d4fd1c 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 = yule_walker_inverse(transfers, noiseCovariance, 50); + + % restricted model iY -> iX - connectivity(iX, iY) = restCovFull(iX) / unCovFull(1, 1) - 1; % source = iY, sink = iX - connectivity(iY, iX) = restCovFull(iY) / unCovFull(2, 2) - 1; % source = iX, sink = iY + % mask for the coefficients of the restricted model + mask = ones(2); + mask(1,2) = 0; - end + % restricted bivariate model (using masked row-by-row solution of + % YW equations) + [transfers_restricted,noiseCovariance_restricted] = yule_walker_mask(R, mask); - % causality in variance - if isfield(inputs, 'lag') && inputs.lag > 0 && any(abs(residual(1, :) - residual(2, :)) > eps) - - % preprocess the residual first - residual = bst_bsxfun(@rdivide, residual, std(residual, [], 2)); - - % bivariate ARCH estimation - [W, C, D, R, S, information, rhoBest] = ... % bivariate ARCH modeling - bst_vgarch('vec', residual, inputs.nTrials, inputs.lag, 0, inputs.flagELM, 3, 'on', 999, [], inputs.rho, false, [], [], [], []); %#ok - - % 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] = 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 @@ -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 = 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 @@ -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_spectral.m b/toolbox/connectivity/bst_granger_spectral.m index 6f2168894..f3605eb7f 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 = yule_walker_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] = 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 - % 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] = 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 @@ -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 = yule_walker_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] = 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 @@ -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/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 +