MVGC demo
Demonstrates typical usage of the MVGC toolbox on generated VAR data for a 5-node network with known causal structure (see var5_test). Estimates a VAR model and calculates time- and frequency-domain pairwise-conditional Granger causalities (also known as the "causal graph"). Also calculates Seth's causal density measure [2].
This script is a good starting point for learning the MVGC approach to Granger-causal estimation and statistical inference. It may serve as a useful template for your own code. The computational approach demonstrated here will make a lot more sense alongside the reference document> [1], which we strongly recommend you consult, particularly Section 3 on design principles of the toolbox. You might also like to refer to the schema of MVGC computational pathways - algorithms A<n> in this demo refer to the algorithm labels listed there - and the Common variable names and data structures section of the Help documentation.
FAQ: Why do the spectral causalities look so smooth?
This is because spectral quantities are calculated from the estimated VAR, rather than sampled directly. This is in accordance with the MVGC design principle that all causal estimates be based on the estimated VAR model for your data, and guarantees that spectral causalities integrate correctly to time-domain causality as theory requires. See [1] for details.
Note: Do not pre-filter your data prior to GC estimation, except possibly to improve stationarity (e.g notch-filtering to eliminate line noise or high-pass filtering to suppress low-frequency transients). Pre-filtering (of stationary data) may seriously degrade Granger-causal inference! If you want (time-domain) GC over a limited frequency range, rather calculate "band-limited" GC; to do this, calculate frequency-domain GCs over the full frequency range, then integrate over the desired frequency band [3]; see smvgc_to_mvgc.
Contents
References
[1] L. Barnett and A. K. Seth, The MVGC Multivariate Granger Causality Toolbox: A New Approach to Granger-causal Inference, J. Neurosci. Methods 223, 2014 [ preprint ].
[2] A. B. Barrett, L. Barnett and A. K. Seth, "Multivariate Granger causality and generalized variance", Phys. Rev. E 81(4), 2010.
[3] L. Barnett and A. K. Seth, "Behaviour of Granger causality under filtering: Theoretical invariance and practical application", J. Neurosci. Methods 201(2), 2011.
Parameters
ntrials = 10; % number of trials nobs = 1000; % number of observations per trial regmode = 'OLS'; % VAR model estimation regression mode ('OLS', 'LWR' or empty for default) icregmode = 'LWR'; % information criteria regression mode ('OLS', 'LWR' or empty for default) morder = 'AIC'; % model order to use ('actual', 'AIC', 'BIC' or supplied numerical value) momax = 20; % maximum model order for model order estimation acmaxlags = 1000; % maximum autocovariance lags (empty for automatic calculation) tstat = ''; % statistical test for MVGC: 'F' for Granger's F-test (default) or 'chi2' for Geweke's chi2 test alpha = 0.05; % significance level for significance test mhtc = 'FDR'; % multiple hypothesis test correction (see routine 'significance') fs = 200; % sample rate (Hz) fres = []; % frequency resolution (empty for automatic calculation) seed = 0; % random seed (0 for unseeded)
Generate VAR test data (A3)
Note: This is where you would read in your own time series data; it should be assigned to the variable X (see below and Common variable names and data structures).
% Seed random number generator. rng_seed(seed); % Get VAR coefficients for 5-node test network. AT = var5_test; nvars = size(AT,1); % number of variables % Residuals covariance matrix. SIGT = eye(nvars); % Generate multi-trial VAR time series data with normally distributed residuals % for specified coefficients and covariance matrix. ptic('\n*** var_to_tsdata... '); X = var_to_tsdata(AT,SIGT,nobs,ntrials); ptoc;
Model order estimation (A2)
% Calculate information criteria up to specified maximum model order. ptic('\n*** tsdata_to_infocrit\n'); [AIC,BIC,moAIC,moBIC] = tsdata_to_infocrit(X,momax,icregmode); ptoc('*** tsdata_to_infocrit took '); % Plot information criteria. figure(1); clf; plot_tsdata([AIC BIC]',{'AIC','BIC'},1/fs); title('Model order estimation'); amo = size(AT,3); % actual model order fprintf('\nbest model order (AIC) = %d\n',moAIC); fprintf('best model order (BIC) = %d\n',moBIC); fprintf('actual model order = %d\n',amo); % Select model order. if strcmpi(morder,'actual') morder = amo; fprintf('\nusing actual model order = %d\n',morder); elseif strcmpi(morder,'AIC') morder = moAIC; fprintf('\nusing AIC best model order = %d\n',morder); elseif strcmpi(morder,'BIC') morder = moBIC; fprintf('\nusing BIC best model order = %d\n',morder); else fprintf('\nusing specified model order = %d\n',morder); end
VAR model estimation (A2)
% Estimate VAR model of selected order from data. ptic('\n*** tsdata_to_var... '); [A,SIG] = tsdata_to_var(X,morder,regmode); ptoc; % Check for failed regression assert(~isbad(A),'VAR estimation failed'); % NOTE: at this point we have a model and are finished with the data! - all % subsequent calculations work from the estimated VAR parameters A and SIG.
Autocovariance calculation (A5)
% The autocovariance sequence drives many Granger causality calculations (see % next section). Now we calculate the autocovariance sequence G according to the % VAR model, to as many lags as it takes to decay to below the numerical % tolerance level, or to acmaxlags lags if specified (i.e. non-empty). ptic('*** var_to_autocov... '); [G,info] = var_to_autocov(A,SIG,acmaxlags); ptoc; % The above routine does a LOT of error checking and issues useful diagnostics. % If there are problems with your data (e.g. non-stationarity, colinearity, % etc.) there's a good chance it'll show up at this point - and the diagnostics % may supply useful information as to what went wrong. It is thus essential to % report and check for errors here. var_info(info,true); % report results (and bail out on error)
Granger causality calculation: time domain (A13)
% Calculate time-domain pairwise-conditional causalities - this just requires % the autocovariance sequence. ptic('*** autocov_to_pwcgc... '); F = autocov_to_pwcgc(G); ptoc; % Check for failed GC calculation assert(~isbad(F,false),'GC calculation failed'); % Significance test using theoretical null distribution, adjusting for multiple % hypotheses. pval = mvgc_pval(F,morder,nobs,ntrials,1,1,nvars-2,tstat); % take careful note of arguments! sig = significance(pval,alpha,mhtc); % Plot time-domain causal graph, p-values and significance. figure(2); clf; subplot(1,3,1); plot_pw(F); title('Pairwise-conditional GC'); subplot(1,3,2); plot_pw(pval); title('p-values'); subplot(1,3,3); plot_pw(sig); title(['Significant at p = ' num2str(alpha)]) % For good measure we calculate Seth's causal density (cd) measure - the mean % pairwise-conditional causality. We don't have a theoretical sampling % distribution for this. cd = mean(F(~isnan(F))); fprintf('\ncausal density = %f\n',cd);
Granger causality calculation: frequency domain (A14)
% Calculate spectral pairwise-conditional causalities at given frequency % resolution - again, this only requires the autocovariance sequence. ptic('\n*** autocov_to_spwcgc... '); f = autocov_to_spwcgc(G,fres); ptoc; % Check for failed spectral GC calculation assert(~isbad(f,false),'spectral GC calculation failed'); % Plot spectral causal graph. figure(3); clf; plot_spw(f,fs);
Granger causality calculation: frequency domain -> time-domain (A15)
% Check that spectral causalities average (integrate) to time-domain % causalities, as they should according to theory. fprintf('\nchecking that frequency-domain GC integrates to time-domain GC... \n'); Fint = smvgc_to_mvgc(f); % integrate spectral MVGCs mad = maxabs(F-Fint); madthreshold = 1e-5; if mad < madthreshold fprintf('maximum absolute difference OK: = %.2e (< %.2e)\n',mad,madthreshold); else fprintf(2,'WARNING: high maximum absolute difference = %e.2 (> %.2e)\n',mad,madthreshold); end