The Neurophysiological Biomarker Toolbox (NBT)

### Part II: VAR models

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));

### Creating a ''dynamics.var'' object

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 $s(t)$. 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 ($B_1, B_2, B_3$), each of dimensions 4 x 4. Property 'VarCoeffs' stores those matrices by concanating them into a big 4×12 matrix $B = [B_1 B_2 B_3]$. So, if you wanted to get the value of matrix $B_2$ 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 $s_j(t)$ towards node $s_i(t)$ requires that the coefficients $b_{ij}(1), ..., b_{ij}(\tau)$ are all zero. Above we just displayed the coefficients for lag $\tau=2$ 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.

### Generating model observations

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.

### Plotting model observations

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.

### What now?

Now you are ready to go to Part III where we will generate fit a VAR model to some real EEG data.