Why donate?

- Tutorials, follow the NBT course

**Follow us **

You are here: Frontpage : Introduction » Tutorials » Granger causality analysis of EEG data » Part II: VAR models

Table of Contents

In this second part of the tutorial we will test whether Granger-causality (GC) indices actually do what they are expected to do. That is, to identify directed information flows between time-series, if the dynamics of the time-series can be (at least approximately) modeled with a Vector Autoregressive Model (VAR).

But first of all you need to run this command to ensure that you start with a perfectly clean MATLAB workspace:

clear classes;

Then you should also run:

import external.eeglab.eegplot;

which will _import_ EEGLAB's function eegplot into your current namespace. The result is that from now on you should be able to run the following command without getting any error:

eegplot(randn(5,1000));

Let's create a model consisting of 4 nodes and order 3 having a pattern of connections between nodes that looks like a path:

varObj = dynamics.var('NbDims', 4, 'Order', 3, 'Topology', 'path');

Now let's examine the properties of our freshly created VAR model:

>> varObj

varObj =

dynamics.var Package: dynamics

Properties: Order: 3 Alpha: 2 VarCoeffs: [4x12 double] Topology: 'path' InnCovCond: [1 1] NbDims: 4 InnCov: [4x4 double] MixingCond: [] MixingMatrix: [] NbMixtures: 4 NoiseCovCond: [] NoiseCov: [] NoiseCovDet: [] Filter: [] InnFilter: [] StreamState: [625x1 uint32] VarSignals: [] Innovations: [] Noise: [] Observations: []

Methods, Superclasses

Recall the formula of a VAR model:

You may already be able to match some of the properties of `varObj`

with
the VAR formula. For instance, property `'Order`

' obviously corresponds to the
model order, i.e. the number of lags in our model (parameter `p`

in the formula).
Also easy is to realize that `'NbDims`

' is just the number of time-series
(or *nodes*) in our model, that is the number of rows of vector . The
other very important property is `'VarCoeffs`

', which logically stores the VAR
model coefficient matrices. Since the model order is `p=3`

we have 3 such
matrices (), each of dimensions `4 x 4`

. Property `'VarCoeffs`

'
stores those matrices by concanating them into a big `4×12`

matrix
. So, if you wanted to get the value of matrix you
could do:

» varObj.VarCoeffs(:,5:8)

ans =

0.1965 0 0.1496 0 0 0.1761 0 0 0 -0.1688 0.0542 0 0.0052 0 0 0.3445

The numbers that you will see are likely to be different as, by default,
a `dynamics.var`

object acquires random coefficient matrices when the
object is created. However, what you will also observe is that there are
at most three off-diagonal entries that are non-zero. Remember that we said
in the lecture that an absence of information flow from node
towards node requires that the coefficients
are all zero. Above we just displayed
the coefficients for lag but if you display the coefficients
for other lags you will find that the same coefficients are zero at
all lags. This means that, in my model, the pattern of information
flow is as depicted below:

Inspect the coefficients of the model that you created and draw a
similar diagram matching the characteristics of your model. Then,
generate a new model with the same order and dimensionality but
having a `loop`

topology. Draw the information flow diagram that
matches this last model as well.

So far we have created a VAR model (i.e. a `dynamics.var`

object), which
means that we have specified the crucial model parameters. However, we
have not yet generated any observations with our model. Let's generate
10000 observations:

>> [s, varObj] = generate(varObj, 10000) s = [1x10000 double] [1x10000 double] [1x10000 double] [1x10000 double] varObj = dynamics.var Package: dynamics Properties: Order: 3 Alpha: 2 VarCoeffs: [4x12 double] Topology: 'path' InnCovCond: [1 1] NbDims: 4 InnCov: [4x4 double] MixingCond: [] MixingMatrix: [] NbMixtures: 4 NoiseCovCond: [] NoiseCov: [] NoiseCovDet: [] Filter: [] InnFilter: [] StreamState: [625x1 uint32] VarSignals: {4x1 cell} Innovations: {4x1 cell} Noise: [] Observations: {4x1 cell}

The variable `s`

that we just created is a
cell array.
A cell array is just a container of objects of (possibly) different types.
This is the fundamental difference with standard MATLAB arrays, as the
latter can only store objects of the same type. You index cell arrays using
the notation `{}`

. For all the examples that we use in this tutorial you
can just convert the output of generate into a standard MATLAB array, which
are easier to work with than cell arrays:

s = cell2mat(s);

Your current workspace should now look something like this:

>> whos Name Size Bytes Class Attributes

s 4x10000 320000 double varObj 1x1 3176 dynamics.var

Note also that when we ran the command:

`[s, varObj] = generate(varObj, 10000)`

we actually updated the _state_ of object `varObj`

. See how properties
`InnCov`

, `VarSignals`

, `Innovations`

and `Observations`

, that
were before empty, now contain some data.

Property `InnCov`

contains the covariance matrix
of the innovation process. The covariance matrix tells how similar the innovation
processes at different network nodes are, i.e. it tells how similar the
components of vector `h(t)`

are. In our case we have:

>> varObj.InnCov

ans =

1.0000 -0.0000 0.0000 0.0000 -0.0000 1.0000 -0.0000 0.0000 0.0000 -0.0000 1.0000 -0.0000 0.0000 0.0000 -0.0000 1.0000

which clearly indicates that the innovation processes at
different network nodes are mutually uncorrelated (the offdiagonal
elements are all zero) and that they all have the same
variance (the diagonal elements are equal). Recall that the
*innovation process* of the VAR model is what we identified in the
lecture slides with the *intrinsic activations* of different
EEG sources, that is with the activity that is truly unique to
each EEG-source and that is not shared with other EEG sources.

Property `Innovations`

contains the innovation process of the model and
property `Observations`

contains the observations. The latter is identical
to the content of variable `s`

in your MATLAB workspace.

Let's see how the observations of our model look like:

eegplot(s, 'srate', 128); % This is equivalent to the command below eegplot(cell2mat(varObj.Observations), 'srate', 128);

The argument `'srate`

' is used to specify the sampling rate of our data.
Of course, the value that I chose (128 Hz) is totally arbitrary as our
data has not been sampled from any analog signal (e.g. from a real EEG) but
they have been artificially generated.

*NOTE:*
As you may notice from inspecting the model observations, our model
generates seemengly random time-series but they do not really look like
EEG. Later in the tutorial we will see how by tuning the model coefficients
we can make our model produce observations that do look like real EEG.