SAILS is a python package for modelling frequency domain power and connectivity in time-series and networks. It provides implementations of a range of autoregressive model fitting, validation and selection routines to describe the linear dependencies in time-series. The spectral content of the fitted models can be explored with within and between channel frequency metrics.

Features

SAILS currently provides:

• Multivariate autoregressive (MVAR) model fits using OLS or Vieira-Morf

• A range of MVAR model diagnostics (Stablility-index, R-square, Durbin-Watson, Percent-Consistency)

• Model order selection via AIC or BIC

• Decomposition of models into oscillatory modes

• Estimation of power spectra using Fourier or Modal transforms

• Wide range of connectivity metrics

  • Transfer Function, Spectral Matrix
  • Coherency, Magnitude Squared Coherence, Phase Coherence
  • Directed Transfer Function and variants
  • Partial Directed Coherence
  • Isolated Effective Coherence

Quick Start

SAILS can be installed from PyPI using pip:

pip install sails

and used to model and describe frequency content in networks of time-series:

import sails
import numpy as np

# Create a simulated signal
sample_rate = 100
siggen = sails.Baccala2001_fig2()
X = siggen.generate_signal(sample_rate=sample_rate,num_samples=1000)

# Fit an autoregressive model with order 3
model = sails.OLSLinearModel.fit_model(X,np.arange(4))

# Compute power spectra and connectivity
freq_vect = np.linspace(0,sample_rate/2)
metrics = sails.FourierMvarMetrics.initialise(model,sample_rate,freq_vect)

Tutorials

Please see our tutorials and example gallery for help getting started with datawoi analysis in SAILS.

Tutorials

The following set of tutorials will lead you through simulating data, fitting models, exploring attributes of the models and performing various plotting routines.

Tutorial 1 - A pink noise system

In this tutorial we will demonstrate how to set up a simple univariate AR model which models a pink noise process. We will use the model to demonstrate how to extract the transfer function using both a Fourier and Modal estimator.

We start by importing the routines we will need from the numpy, sails and matplotlib. We set the ggplot style for all plots in the tutorials.

import numpy as np
from sails import generate_pink_roots, root_plot, model_from_roots
import matplotlib.pyplot as plt

plt.style.use('ggplot')

We now define the meta-parameters of our system. We will arbitrarily pretend that we are using a sampling rate of 100Hz. From this, we can compute our Nyquist frequency, and a vector of frequencies at which we will evaluate measures such as our Fourier estimator. The variable freq_vect will come up often during our tutorials. In this case, we will estimate 64 frequencies which are linearly spaced between 0Hz and Nyquist. You can alter freq_vect and examine how this affects the results of the later estimations.

sample_rate = 100
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

As mentioned above, for our first system we will emulate a system which behaves as a pink noise process. Whereas we would normally learn our model from some data, for this example we will construct a model from the polynomial form of the AR model. A function (generate_pink_roots()) has been provided which calculates a set of roots which will generate the appropriate roots:

roots = generate_pink_roots(1)

Given the roots of the model, we can plot these in a Z-plane representation to examine them. A function (root_plot()) is provided to make this straightforward.

ax = root_plot(roots)
ax.figure.show()

The plot shows frequency increasing counterclockwise from the x-axis. Nyquist frequency is on the negative x-axis. For real systems, the structure of the pole plot will be mirrored across the x-axis.

As mentioned above, we would normally set up a model by learning the parameters from data. During this exploratory tutorial, we can analytically create the model from the roots. The model_from_roots() function will perform this task for you.

m = model_from_roots(roots)

No matter how we create our model, it will be a subclass of AbstractLinearModel. All subclasses of AbstractLinearModel provide a number of guarantees so that they can all be used with the various estimator routines we discuss below. We will discuss fitting a model from data in a later tutorial.

Now that we have our model, we can extract the transfer function (H) using two different methods. The first method extracts the transfer function using a Fourier transform of the parameters.

from sails import FourierMvarMetrics

F = FourierMvarMetrics.initialise(m, sample_rate, freq_vect)
F_H = F.H

The second method using a modal decomposition of the parameter matrix to break the parameters apart into, hopefully, interpretable components. We first of all create a ModalMvarMetrics and then use the MvarModalDecomposition object inside of this to extract the transfer function for each mode. Each mode will consist of either a pair of poles in the complex plane or a single real mode. The per_mode_transfer_function() will take this into account for you.

from sails import ModalMvarMetrics

M = ModalMvarMetrics.initialise(m, sample_rate, freq_vect)
M_H = M.modes.per_mode_transfer_function(sample_rate, freq_vect)

We can plot our modes in both forms to examine the relationship between them. Firstly, we sum the modal transfer function across all modes to recover the Fourier based transfer function (to within computational precision).

fourier_H = np.abs(F_H).squeeze()
modal_H = np.abs(M_H.sum(axis=1)).squeeze()

# Check the two forms are equivalent
equiv = np.allclose( fourier_H, modal_H )
ssdiff = np.sum( (fourier_H - modal_H)**2 )
print('Fourier and Modal transfer functions equivalent: {0}'.format(equiv))
print('Sum-square residual difference: {0}'.format(ssdiff))

# Plot our fourier and modal spectra

f2 = plt.figure()

plt.plot(freq_vect, fourier_H, 'o');
plt.plot(freq_vect, modal_H);

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

plt.legend(['Fourier H', 'Modal H'])

plt.savefig('tutorial1_2.png', dpi=300)
f2.show()

Finally, we can see how each mode contributes to the overall transfer function shape by plotting each mode separately. Each mode contributes a single uni-modal resonance to the spectrum. In this case there are no clear peaks in the spectrum, just the 1/f shape - as such, all the mode peaks sum together to make a smooth spectrum. Note that we’re plotting the modes on a log y-scale as some modes make a very small contribution to the overall transfer function.

fourier_H = np.abs(F_H).squeeze()
modal_H = np.abs(M_H).squeeze()

# Plot our fourier and modal spectra

f2 = plt.figure()

plt.semilogy(freq_vect, fourier_H, 'o');
plt.semilogy(freq_vect, modal_H);

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

legend = ['Mode' + str(ii) for ii in range(modal_H.shape[1])]
plt.legend(['Fourier H'] + legend)

plt.savefig('tutorial1_3.png', dpi=300)
f2.show()

Tutorial 2 - A more interesting oscillation

In this tutorial, we will extend our example to include oscillations on top of the background pink noise. We begin in the same way as with the previous tutorial.

import numpy as np
from sails import generate_pink_roots, root_plot, model_from_roots
import matplotlib.pyplot as plt

plt.style.use('ggplot')

sample_rate = 100
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

roots = generate_pink_roots(1)

We now add in energy at two frequencies and plot the poles.

roots[3:5] = roots[3:5]*1.1

roots[7:9] = roots[7:9]*1.1

ax = root_plot(roots)

ax.figure.show()

In this plot we can see that in addition to our pink noise mode on the x-axis (compare back to the equivalent plot in the first tutorial), we have two poles which now sit closer to the unit circle - these correspond to the additional energy which we have placed in the signal.

We can now continue to set up our model from our roots, generate the Fourier and Modal transfer functions and plot them as before:

from sails import FourierMvarMetrics, ModalMvarMetrics

m = model_from_roots(roots)

F = FourierMvarMetrics.initialise(m, sample_rate, freq_vect)
F_H = F.H

M = ModalMvarMetrics.initialise(m, sample_rate, freq_vect)
M_H = M.modes.transfer_function(sample_rate, freq_vect, sum_modes=True)

# Plot our fourier and modal spectra
f2 = plt.figure()

plt.plot(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.plot(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f2.show()

In comparison to our plot in the first tutorial, there are (as expected) two peaks in the Fourier transfer function at around 8Hz and 17Hz. We can also see that each of these two peaks has been extracted into a separate mode in the modal decomposition. This can be seen slightly more clearly by using a log scale for the y axis:

f3 = plt.figure()

plt.semilogy(freq_vect, np.abs(F_H).squeeze(), 'o');

plt.semilogy(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')

plt.ylabel('Frequency Response')

f3.show()

We can also go on to extract the magnitude of the eigenvalues for each of the poles:

ev = np.abs(M.modes.evals)

Note that we have an eigenvalue for each of the poles, not each of the modes. As the eigenvalue will be identical for both of the poles in a given mode (where the mode consists of a pole pair), we can reduce down to examining one pole for each mode. To get the relevant indices, we can extract the information from the mode_indices property on the MvarModalDecomposition object. We can then extract other interesting pieces of information such as the period of the oscillation in the same way:

idx = [i[0] for i in M.modes.mode_indices]

print(ev[idx])

[[ 0.96393242]
 [ 0.83308007]
 [ 0.88077379]
 [ 0.78215297]
 [ 0.84650308]
 [ 0.7603617 ]
 [ 0.75345934]
 [ 0.74824023]
 [ 0.74435422]
 [ 0.73893252]
 [ 0.73980532]
 [ 0.74158676]]

print(M.modes.peak_frequency[idx])

[[  0.        ]
 [  4.37353106]
 [  8.83055516]
 [ 13.2106547 ]
 [ 17.56247546]
 [ 21.90065104]
 [ 26.23121568]
 [ 30.55715939]
 [ 34.88016321]
 [ 47.84046803]
 [ 43.52119577]
 [ 39.20127164]]

From this, we can see that the modes which fit our two oscillations are found at 8.8Hz and 17.6Hz.

Tutorial 3 - Fitting real univariate data

In the previous two tutorials we set up our system using the polynomial representation. In most cases, we will want to learn the parameters of our model from real data. In this section, we cover how this is done.

We will be using some data from a single MEG Virtual Electrode reconstruction. The data has been down-sampled from the original sampling rate for the purposes of demonstration and ease of analysis.

We start with some imports and finding the location of our example data:

from os.path import join
import numpy as np
import matplotlib.pyplot as plt

from sails import find_example_path, root_plot, model_from_roots

plt.style.use('ggplot')

ex_dir = find_example_path()

The original MEG data was sampled at 2034.51Hz and has been down-sampled by 24 times giving us a sampling rate of just under 85Hz. The data has also been Z-scored. The data is stored in an HDF5 file in the example data directory as meg_single.hdf5; the data is stored as the dataset X.

import h5py

sample_rate = 2034.51 / 24

nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

X = h5py.File(join(ex_dir, 'meg_single.hdf5'), 'r')['X'][...]

print(X.shape)

(1, 30157, 1)

The form of the data for the fitting routines is (nsignals, nsamples, ntrials). In the current case, we are performing a univariate analysis so we only have one signal. We have just over 30,000 data points and one one trial.

We are now in a position to set up our model. We will use the Vieira-Morf algorithm to fit our model, and fit a model of order 19 using the delay_vect argument. We will discuss model order selection in a later tutorial.

from sails import VieiraMorfLinearModel

delay_vect = np.arange(20)

m = VieiraMorfLinearModel.fit_model(X, delay_vect)

print(m.order)

19

We can now compute some diagnostics in order to check that our model looks sensible. We can start by looking at the R-squared value:

diag = m.compute_diagnostics(X)

print(diag.R_square)

0.41148872414

We can see that we are explaining just over 40% of the variance with our model which, given that we are modelling human MEG data collected over roughly 6 minutes, is reasonable.

The diagnostics class also gives us access to various other measures:

print(diag.AIC)

-25563.4038876

print(diag.BIC)

-25396.883104

print(diag.LL)

-25603.4038876

print(diag.DW)

[ 2.00103186]

print(diag.SI)

0.976853801393

print(diag.SR)

0.0663411856574

In turn, these are the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Log Likelihood (LL), Durbin-Watson coefficient (DW), Stability Index (SI) and Stability Ratio (SR). It is also possible to access the Percent Consistency (PC), although this is not computed by default due to it being memory intensive - you can compute this using the percent_consistency() routine.

As in our previous examples, we can extract our metrics and plot the transfer functions using both the Fourier and Modal methods as we have previously done:

from sails import FourierMvarMetrics, ModalMvarMetrics

F = FourierMvarMetrics.initialise(m, sample_rate, freq_vect)
F_H = F.H

M = ModalMvarMetrics.initialise(m, sample_rate, freq_vect)
M_H = M.modes.per_mode_transfer_function(sample_rate, freq_vect)

# Plot our fourier and modal spectra

f2 = plt.figure()

plt.plot(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.plot(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f2.show()

f3 = plt.figure()

plt.semilogy(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.semilogy(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f3.show()

In our previous examples, the model was defined by the structure of the polynomial, and we could analytically write down the form of the poles. In this example, we have learned the parameters from data. We may want to look at the plot of the system roots:

ax = M.modes.pole_plot()
ax.figure.show()

As previously, we can also go on to extract the magnitude of the eigenvalues and period of each of the modes:

ev = np.abs(M.modes.evals)

idx = [i[0] for i in M.modes.mode_indices]

print(ev[idx])

[[ 0.9768538 ]
 [ 0.88790649]
 [ 0.83392372]
 [ 0.82257801]
 [ 0.7869429 ]
 [ 0.79447331]
 [ 0.80151847]
 [ 0.81061091]
 [ 0.81709057]
 [ 0.80730369]]

print(M.modes.peak_frequency[idx])

[[  0.        ]
 [  8.88121272]
 [  4.51372792]
 [ 13.89286835]
 [ 18.15521077]
 [ 22.01471814]
 [ 26.47871083]
 [ 30.9324283 ]
 [ 35.59229698]
 [ 40.1740829 ]]

From this, we can see that the mode which primarily fits this data is an alpha mode at 8.9Hz.

Tutorial 4 - Exploring model order

In this tutorial, we will examine how to look at different model orders in order to determine the most appropriate model.

The first section of this tutorial should be familiar to you from previous tutorials. We start by loading in some data:

from os.path import join

import h5py

import numpy as np
import matplotlib.pyplot as plt

from sails import find_example_path, VieiraMorfLinearModel

plt.style.use('ggplot')

ex_dir = find_example_path()

sample_rate = 2034.51 / 24
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

X = h5py.File(join(ex_dir, 'meg_single.hdf5'), 'r')['X'][...]

We will use the AIC values to examine how appropriate our model is. We therefore set up an empty list to store our AIC values in and fit a model for a range of different orders. We can then plot up the values of AIC.

AICs = []

orders = []

for delays in range(2, 35):
    delay_vect = np.arange(delays)
    m = VieiraMorfLinearModel.fit_model(X, delay_vect)
    diag = m.compute_diagnostics(X)
    orders.append(m.order)
    AICs.append(diag.AIC)

f = plt.figure()

plt.plot(orders, AICs, 'o');

plt.xlabel('Model Order')
plt.ylabel('AIC')

f.show()

The same principle can be used for any of the other measures such as BIC and R-squared. You should always be aware that the change the number of parameters in the model for different orders means that certain measures (such as R-squared) will always increase for increasing model order.

R_squares = []

orders = []

for delays in range(2, 35):
    delay_vect = np.arange(delays)
    m = VieiraMorfLinearModel.fit_model(X, delay_vect)
    diag = m.compute_diagnostics(X)
    orders.append(m.order)
    R_squares.append(diag.R_square)

f2 = plt.figure()

plt.plot(orders, R_squares, 'o');

plt.xlabel('Model Order')
plt.ylabel('$R^2$')

f2.show()

Tutorial 5 - MVAR Connectivity Estimation

In this tutorial, we will explore a range of connectivity estimators in a simulated network.

We start by importing sails and defining some meta-parameters as we did in previous tutorials.

import numpy as np

import sails

sample_rate = 128

We will use an example network from Baccala & Sameshima 2001. This is defined within the simulation module in Sails. We can import a signal generator based on figure 2 from Baccala & Sameshima 2001.

siggen = sails.simulate.Baccala2001_fig2()

The siggen object contains the autoregressive parameters defining the network and can create a model containing these parameters which we can use in further analyses. Here we create a model containing the ‘true’ autoregressive parameters and use it to compute the connectivity metrics.

m = siggen.generate_model()

freq_vect = np.linspace(0, sample_rate/2, 36)

F = sails.FourierMvarMetrics.initialise(m, sample_rate, freq_vect)

F is an object containing methods to compute a range of frequency domain connectivity metrics. Each metric is evaluated across the range of frequencies_defined in freq_vect and can be plotted using the plot_vector function in sails. Next we plot the Magnitude Squared Coherence estimated from our simulation parameter.

fig = sails.plotting.plot_vector(np.abs(F.S), freq_vect, diag=True)

fig.show()

This will generate a matrix of plots, each plot represents the coherence as a function of frequency between the node specified in the column and row labels. In this case, we find that all our nodes are strongly coherence with each other. This is as the coherence does not distinguish between direct and indirect connections. For example, nodes 1 and 5 are only connected through node 4 yet the coherence still shows a connection. The coherence is also symmetric, that is the connection from node 1->2 is the same as node 2->1.

Next we plot the Directed Transfer Function which is a directed measure that is able to show when connections are not symmetrical

fig = sails.plotting.plot_vector(F.directed_transfer_function, freq_vect, diag=True)

fig.show()

The Directed Transfer Function shows far fewer connections than the Magnitude Squared Coherence. We can now see that the connections between node 1 and the rest of the nodes are asymmetrical. This means that node 1 is driving the others. The interaction between nodes 4 and 5 is now also isolated. A remaining issue is that the Directed Transfer Function is still sensitive to indirect connections, as we can see by the power in the subplot between node 1 and 5. The Partial Directed Coherence aims to address this problem.

fig = sails.plotting.plot_vector(F.partial_directed_coherence, freq_vect, diag=True)

fig.show()

The Partial Directed Coherence now shows only the direct connections within our network. We retain our frequency resolution and the sensitivity to asymmetrical connections. There are many other MVAR derived connectivity metrics available within sails with different properties and sensitivities, these include:

• Coherency (sails.mvar_metrics.coherency(), sails.mvar_metrics.AbstractMVARMetrics.coherency)
• Imaginary Coherence (sails.mvar_metrics.AbstractMVARMetrics.imaginary_coherence)
• Phase Coherence (sails.mvar_metrics.AbstractMVARMetrics.phase_coherence)
• Magnitude Squared Coherence (sails.mvar_metrics.AbstractMVARMetrics.magnitude_squared_coherence)
• Partial Coherence (sails.mvar_metrics.partial_coherence(), sails.mvar_metrics.AbstractMVARMetrics.partial_coherence)
• Directed Transfer Function (sails.mvar_metrics.directed_transfer_function(), sails.mvar_metrics.AbstractMVARMetrics.directed_transfer_function)
• Full Frequency Directed Transfer Function (sails.mvar_metrics.AbstractMVARMetrics.ff_directed_transfer_function)
• Directed Directed Transfer Function (sails.mvar_metrics.AbstractMVARMetrics.d_directed_transfer_function)
• Partial Directed Coherence (sails.mvar_metrics.partial_directed_coherence(), sails.mvar_metrics.AbstractMVARMetrics.partial_directed_coherence)
• Isolated Effective Coherence (sails.mvar_metrics.isolated_effective_coherence(), sails.mvar_metrics.AbstractMVARMetrics.isolated_effective_coherence)

In the second part of this tutorial we will look at fitting and MVAR model and the Partial Directed Coherence to simulated data, rather than from the ‘true’ model.

We can generate data from our simulated model using the sails.simulate:AbstractSigGen.generate_signal() method and specifying the sample_rate and number of samples to generate

X = siggen.generate_signal(sample_rate=128, num_samples=640)

X is a (nchannels x nsamples) array containing our simulated data. We can plot X using matplotlib

import matplotlib.pyplot as plt

plt.figure()

for ii in range(5):
    plt.plot(X[ii, :] + (ii * 10))

plt.show()

We now have a figure containing 5 time-series from our simulation. We can see there is an oscillation by eye and that some of the time-series vary together more than others.

We can fit a model to the simulated data and compute connectivity metrics as we did in previous tutorials.

delay_vect = np.arange(4)

m = sails.VieiraMorfLinearModel.fit_model(X, delay_vect)

F = sails.FourierMvarMetrics.initialise(m, sample_rate, freq_vect)

diag = m.compute_diagnostics(X)

We check that our model is fitting well by interrogating the diagnostics. Here we see that we are explaining around 56% of the total variance in the signal and that our model is stable (diag.SI = .91).

Let’s compare the Partial Directed Coherence from our fitted model to the Partial Directed Coherence from the ‘true’ model.

m0 = siggen.generate_model() # This is our true model

F0 = sails.FourierMvarMetrics.initialise(m0, sample_rate, freq_vect)

pdc = np.concatenate((F0.partial_directed_coherence, F.partial_directed_coherence), axis=3)

fig = sails.plotting.plot_vector(pdc, freq_vect, diag=True, line_labels=['True', 'Fitted'])

fig.show()

The resulting figure shows the nodes by nodes matrix of subplots containing the PDC estimates. We can see that our model is doing a pretty good job approximating the true pattern of connectivity. There may be some false-positive connections which show power for the fitted model but not for the true model.

Try re-running the simulation with a higher or lower number of samples in the time series. You should see that the estimation starts to really break down (lots of false positives and a distorted spectrum shape) when we have too few samples (e.g. num_samples = 128) and becomes nearly perfect when we have a very long time-series (e.g. num_sample = 2048)

Tutorial 6 - Sliding window univariate model estimation

In this tutorial, we will examine fitting multiple models using a sliding window approach.

For this tutorial, we will use the same MEG example data which we have used in previous tutorials.

We start by importing our modules and finding and loading the example data.

import os
from os.path import join

import h5py

import numpy as np
import matplotlib.pyplot as plt

import sails

plt.style.use('ggplot')

SAILS will automatically detect the example data if downloaded into your home directory. If you’ve used a different location, you can specify this in an environment variable named SAILS_EXAMPLE_DATA.

# Specify environment variable with example data location
# This is only necessary if you have not checked out the
# example data into $HOME/sails-example-data
os.environ['SAILS_EXAMPLE_DATA'] = '/path/to/sails-example-data'

# Locate and validate example data directory
example_path = sails.find_example_path()

# Load data using h5py
data_path = os.path.join(sails.find_example_path(), 'meg_occipital_ve.hdf5')
X = h5py.File(data_path, 'r')['X'][:, 122500:130000, :]

X = sails.utils.fast_resample(X, 5/2)

sample_rate = 250 / (5/2)
nyq = sample_rate / 2.

time_vect = np.linspace(0, X.shape[1] // sample_rate, X.shape[1])
freq_vect = np.linspace(0, nyq, 64)

As a reminder, our data is (nsignals, nsamples, ntrials):

print(X.shape)

(1, 3000, 1)

The idea behind a sliding window analysis is to fit a separate MVAR model to short windows of the data. For the purposes of this example, we will fit consecutive windows which are 200 samples long. We can do this by simply extracting the relevant data and fitting the model in the same way that we have done so before. We are using 10 delays in each of our models. After computing the model, we can then calculate some of our diagnostics and examine, for example, the \(R^2\) value. We could also choose to examine our Fourier MVAR metrics.

from sails import VieiraMorfLinearModel, FourierMvarMetrics

delay_vect = np.arange(10)

m1 = VieiraMorfLinearModel.fit_model(X[:, 0:200, :], delay_vect)

d1 = m1.compute_diagnostics(X[:, 0:200, :])

print(d1.R_square)

0.6695709340724934

We can then repeat this procedure for the next, overlapping, window:

delay_vect = np.arange(10)

m2 = VieiraMorfLinearModel.fit_model(X[:, 1:201, :], delay_vect)

d2 = m1.compute_diagnostics(X[:, 1:201, :])

print(d2.R_square)

0.6567768337622513

This is obviously a time-consuming task, so we provide a helper function which will take a set of time-series data and compute a model and diagnostic information for each window. The sails.modelfit.sliding_window_fit() function takes at least five parameters:

1. the class to use for fitting (here, we use sails.modelfit.VieiraMorfLinearModel)
2. the data to which the models are being fitted
3. the delay vector to use for each model
4. the length of each window in samples
5. the step between windows. For example, if this is set to 1 sample, consecutive overlapping windows will be used. If this is set to 10, each consecutive window will start 10 samples apart. If this is equal to the length of the window or greater, the windows will not overlap.

The sliding_window_fit routine returns specially constructed model in which the parameters matrix is set up as (nsources, nsources, ndelays, nwindows). In other words, there is a model (indexed on the final dimension) for each of the windows. A ModelDiagnostics class is also returned which contains the diagnostic information for all of the windows. Use of the sliding_window_fit() routine is straightforward:

from sails import sliding_window_fit

M, D = sliding_window_fit(VieiraMorfLinearModel, X, delay_vect, 100, 1)

model_time_vect = time_vect[M.time_vect.astype(int)]

Once we have constructed our models, we can go ahead and construct our FourierMvarMetrics class which will allow us to, for example, extract the transfer function for each window:

F = FourierMvarMetrics.initialise(M, sample_rate, freq_vect)

We will now produce a composite plot which illustrates how the behaviour of the system evolves over time. We are going to limit ourselves to the first 10000 data points and windows.

f1, axes = plt.subplots(nrows=5, ncols=1, figsize=(12, 8))

On the top row, we plot our data.

plt.subplot(5, 1, 1)

plt.plot(time_vect, X[0, :, 0])

plt.xlim(0, 30)
plt.grid(True)

plt.ylabel('Amplitude')

On the next two rows, we plot our transfer function for the first 10000 windows:

plt.subplot(5, 1, (2, 3))

plt.contourf(model_time_vect, F.freq_vect[3:], np.sqrt(np.abs(F.PSD[0, 0, 3:, :])))

plt.xlim(0, 30)
plt.grid(True)

plt.ylabel('Frequency (Hz)')

Underneath the transfer function, we examine the stability index (in red) and the \(R^2\) value (in blue):

plt.subplot(5, 1, 4)

plt.plot(model_time_vect, D.SI)

plt.plot(model_time_vect, D.R_square)

plt.xlim(0, 30)
plt.grid(True)

plt.ylabel('SI / $R^2$')

On the bottom row, we plot the AIC value for each of the models:

plt.subplot(5, 1, 5)

plt.plot(model_time_vect, D.AIC)

plt.xlim(0, 30)
plt.grid(True)

plt.xlabel('Time (samples)')

plt.ylabel('AIC')

Finally, we can look at our overall figure:

f1.tight_layout()

f1.show()

Tutorial 7 - Plotting Helpers

In this tutorial, we will explore the plotting helper functions which are available for use in sails; primarily for plotting netmats and circular connectivity pots.

For this tutorial, we will use some example plot definition files which are provided with sails.

We start by importing our modules and finding and finding our example path

from os.path import join

from sails import find_support_path

support_dir = find_support_path()

group_csv = join(support_dir, 'aal_cortical_plotting_groups.csv')
region_csv = join(support_dir, 'aal_cortical_plotting_regions.csv')

The two files which are imported above describe the layout of the plots which we will create. The example files given are for a specific cortical subset of the AAL atlas which we have used in some of our own work; this parcellation contains 78 regions; 39 per hemisphere.

The groups.csv file contains a list of the groups of regions which we will plot. In this case, we have divided the 78 regions into 10 groups:

• Frontal left / right
• Medial left / right
• Temporal left / right
• Parietal left / right
• Occipital left / right

The regions.csv file describes each ROI and places it within these groups. It also provides the indexing information so that we can match between our netmats and the plots.

More details of both of these files can be found in the docstring for the from_csv_files function of the sails.CircosHandler class.

We will now build a connectivity diagram using this structure.

Setting up the handler

We start by setting up a handler which parses the CSV files:

from sails import CircosHandler

c = CircosHandler.from_csv_files(group_csv, region_csv)

There is one additional argument available to the from_csv_files function; this argument is analysis_column. By default, the order in which our analysed data is stored comes from the AnalysisID column of the regions.csv spreadsheet. You can alter which column is used for this by setting the analysis_column argument to the function. This is useful in cases where you have data in different orders from different analyses but where you want to produce plots with the same ordering with little effort.

Extracting Region Ordering and Indices

We can use use this handler to extract various orderings of our ROIs and the indices needed to re-order data into them.

# Show the order in circular form (clockwise) of the groups
print(c.circular_groups)

['frontal_L', 'medial_L', 'temporal_L', 'parietal_L', 'occipital_L', 'occipital_R', 'parietal_R', 'temporal_R', 'medial_R', 'frontal_R']

# Show the order in netmat form (top to bottom) of the groups
print(c.netmat_groups)

['frontal_L', 'medial_L', 'temporal_L', 'parietal_L', 'occipital_L', 'frontal_R', 'medial_R', 'temporal_R', 'parietal_R', 'occipital_R']

# Show the last 10 regions in circular order and the corresponding
# indices into our data
print(c.circular_regions[-10:])
print(c.circular_indices[-10:])

['Frontal_Sup_Medial_R', 'Frontal_Inf_Tri_R', 'Frontal_Inf_Oper_R', 'Frontal_Mid_R', 'Frontal_Sup_R', 'Frontal_Inf_Orb_R', 'Frontal_Mid_Orb_R', 'Frontal_Med_Orb_R', 'Frontal_Sup_Orb_R', 'Rectus_R']
[28, 19, 15, 25, 31, 17, 24, 21, 30, 61]

# Show the last 10 regions in netmat order and the corresponding
# indices into our data
print(c.netmat_regions[-10:])
print(c.netmat_indices[-10:])

['SupraMarginal_R', 'Rolandic_Oper_R', 'Precuneus_R', 'Occipital_Sup_R', 'Occipital_Mid_R', 'Occipital_Inf_R', 'Calcarine_R', 'Cuneus_R', 'Lingual_R', 'Fusiform_R']
[67, 63, 59, 45, 43, 41, 5, 13, 39, 33]

Circular Plots

To generate circular plots, we use the program circos. You will need to have the program installed on your computer in order to generate such plots. circos is available from http://circos.ca/software/. On Debian and similar, you can simply apt install circos circos-data libsvg-perl.

Note that if you are creating plots for publication using Circos, you should cite the relevant publication: Krzywinski, M. et al. Circos: an Information Aesthetic for Comparative Genomics. Genome Res (2009) 19:1639-1645.

Generating a karyotype

Circos uses the term karyotype to describe the ordering of chromosomes and bands within them. In our case, we are using this to describe groups of regions and individual ROIs respectively. We only need to generate a karyotype file once for each layout; we do not need a different karyotype file for each individual plot.

To generate a karyotype file, we open a file and write the contents of the karyotype() function into it.

f = open('aal_karyotype.txt', 'w')
f.write(c.karyotype())
f.close()

In our example case, the start of the file will look like this:

chr - frontal_L 1 0 12000 red
chr - medial_L 2 0 6000 yellow
chr - temporal_L 3 0 6000 purple
chr - parietal_L 4 0 8000 green
chr - occipital_L 5 0 7000 blue
chr - occipital_R 6 0 7000 vvdblue
chr - parietal_R 7 0 8000 vvdgreen
chr - temporal_R 8 0 6000 vvdpurple
chr - medial_R 9 0 6000 vvdyellow
chr - frontal_R 10 0 12000 vvdred
band frontal_L Rectus_L Rectus_L 0 1000 grey
band frontal_L Frontal_Sup_Orb_L Frontal_Sup_Orb_L 1000 2000 grey

We will use this file in the next section to make our plots.

Generating a connectivity plot based on the karyotype

We start by generating some controlled data with connections between only a few regions

import numpy as np

data = np.zeros((78, 78)) * np.nan

# Add a strong positive connection between Amygdala_L and Insula_R
data[0, 37] = 10.0
# Add a weaker negative connection between Cuneus_R and Heschl_R
data[13, 35] = -4.0

Note that when the data values are used in the circular plotting routines, they will be used as line widths in pixels.

We now generate a set of circos configuration files from this data. We also need to pass it our karyotype file name and the output base name.

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out')

The code above assumes that you have the circos config files installed in /etc/circos. If you have them in another location, pass the circos_path variable to the gen_circos_files routine, e.g.:

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out',
                   circos_path='/usr/local/etc/circos')

This will generate two files: aal_test_out.conf and aal_test_out.txt. The former is the circos configuration file and the latter is the file which contains the actual information about the connections. In our case, we can see that the latter only contains two lines; one for each of our two connections.

Manually Generating the plot

To manually generate the plot from the configuration files, we use a normal shell and run the circos command.

circos -conf aal_test_out.conf

You may get an error which contains the following:

*** CIRCOS ERROR ***

cwd: /tmp

command: /usr/bin/circos -conf aal_test_out.conf

CONFIGURATION FILE ERROR

This is because of a problem with circos finding some of its configuration files. You can fix this by running the following commands from your shell (whilst in the directory containing the files). This assumes that the circos config files are in /etc/circos:

mkdir etc
    ln -s /etc/circos/tracks etc/

You should then re-run the circos command.

Two files will be created: test_out.png and test_out.svg. The image should look like this:

Automatically Generating the plot

Instead of going to the effort of the above, there is a helper routine gen_circos_plot which will create a temporary directory, generate all of the configuration files, run circos and then copy the plots to your final destination. It can be used directly in place of gen_circos_files. If the circos call fails, an exception will be raised which will include the contents of stdout and stderr so that you can investigate the problem.

Note that you do not pass a karyotype file name to the gen_circos_plot routine as it will generate a karyotype file in the temporary directory for you.

c.gen_circos_plot(data, 'aal_test_out',
                  circos_path='/usr/local/etc/circos')

Modifying link colours

By default, links with a negative value will be shown in blue, links with a positive value in red and links with a strictly zero value in black. You can modify this in two ways:

1. By setting a tuple for (zero_colour, negative_colour, positive_colour)
2. By passing a numpy matrix with dtype object containing a string for each connection. In this case, you can also pass a defaultcolor which will be used if a matrix entry is None

An example of using the tuple syntax:

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out_2',
                   linkcolors=('black', 'green', 'purple'))

and an example of using the matrix syntax:

# As an example here, we set the colour of the Amygdala/Insula link
# explicitly and set the other link using the default color syntax
colors = np.empty((78, 78), dtype=object)

colors[0, 37] = 'yellow'

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out_3',
                   linkcolors=colors,
                   defaultlinkcolor='orange')

Netmat Plots

In this section, we will discuss how to use the netmat plotting routines which come as part of the CircosHandler above. Further down there is documentation on how to use the raw plot_netmat routine.

We assume that we are in the same session as above and that c still represents a CircosHandler object.

# Set up some netmat data
netdata = np.zeros((78, 78))

# Link from Calcarine_L (4) to Heschl_R (35) with a positive score
netdata[4, 35] = 1.0

# Link from Calcarine_L (4) to Heschl_L (34) with a negative score
netdata[4, 34] = -1.0

# Link from Frontal_Mid_L (22) to Fusiform_L (32) with a double positive socre
netdata[22, 32] = 2.0

# Create our plot - a figure and axes will be created as we have
# not supplied any

# Note that, as demonstrated here, you can use named arguments to pass
# extra options into the pcolormesh call.  See the docstring of
# the CircosHandler.plot_netmat function for more help
ax = c.plot_netmat(netdata, label_fontsize=4, cmap='RdBu_r', vmin=-2.0, vmax=2.0)

Raw Netmat Plots

If you do not wish to set up a full set of CSV files etc to plot your netmats, you can call the sails.plotting.plot_netmat routine having prepared your arguments yourself.

You will need to set up colour mappings and lists as well as (optionally) labels for your regions. An example is shown below which uses the same data as above but changes the colour layout slightly:

from sails import plot_netmat

    # These original colours were taken as RGB tuples from circos
orig_colors = {'red':        (217, 120, 99),
               'orange':     (217, 144, 89),
               'yellow':     (220, 213, 103),
               'purple':     (155, 152, 183),
               'green':      (127, 180, 128),
               'blue':       (121, 166, 193),
               'vvdred':     (152, 49, 58),
               'vvdorange':  (143, 79, 52),
               'vvdyellow':  (178, 170, 49),
               'vvdpurple':  (99, 62, 139),
               'vvdgreen':   (49, 109, 82),
               'vvdblue':    (54, 95, 148)}

# We need to convert them to to matplotlib compatible RGB 0-1 tuples
colors = {}

    for cname, cval in orig_colors.items():
        colors[cname] = [(x / 255.0) for x in cval]

    # This is the layout of our AAL parcellation division
    cnames = ['red']       * 12 + \
         ['yellow']    * 6 + \
         ['purple']    * 6 + \
         ['green']     * 8 + \
         ['blue']      * 7 + \
         ['vvdred']    * 12 + \
         ['vvdyellow'] * 6 + \
         ['vvdpurple'] * 6 + \
         ['vvdgreen']  * 8 + \
         ['vvdblue']   * 7

tick_pos = [12, 14, 18, 24, 32, 39, 51, 53, 57, 63, 71]

    # Some temporary labels for testing
    labels = ['R{}'.format(x) for x in range(1, 79)]

    # Note that in reality, we would have to re-order our data to be
# in the correct order for the netmat here (in the above example,
# this is handled by the CircosHandler class).
# Here we leave it unordered, which is why the plot ends up
# different (and wrong given what we claimed we were doing with
# the data above)

    plot_netmat(netdata, colors=colors, cnames=cnames, tick_pos=tick_pos,
        labels=labels, label_fontsize=4, cmap='RdBu_r', vmin=-2.0, vmax=2.0)

Tutorial 8 - Using A Custom Model Fit

SAILS provides implementations of several algorithms for fitting autoregressive models but it is straightforward to create a custom class which implements a new model fit or uses one from another package.

This tutorial will outline how to create a custom model fit class using the Vector Autoregression class from statsmodels.tsa. We start by importing SAILS and creating a simulated time series to model.

import sails

# Create a simulated signal
sample_rate = 100
siggen = sails.Baccala2001_fig2()
X = siggen.generate_signal(sample_rate=sample_rate,num_samples=1000)

We then fit a model using Ordinary Least Squared as implemented in SAILS.

# Fit an autoregressive model with order 3
sails_model = sails.OLSLinearModel.fit_model(X, np.arange(4))

Before we continue, we should note that you can pass an existing sklearn estimator (for example, sklearn.linear_model.LinearRegression as the estimator parameter to the fit_model function of the OLSLinearModel class. If you do this, you should not fit the intercept in the model. For instance:

import sklearn.linear_model

# Fit an autoregressive model using SKLearn's LinearRegressor
estimator = sklearn.linear_model.LinearRegression(fit_intercept=False)
sails_model = sails.OLSLinearModel.fit_model(X, np.arange(4), estimator=estimator)

The above will give the same answers as the default method (which calculates the parameters using the normal equations). You can, however, extend this approach to use, for instance, ridge or lasso-regression using the relevant classes (sklearn.linear_model.Ridge or sklearn.linear_model.Lasso).

To go beyond what is available using the previous options, we can create a new model fit class based on sails.modelfit.AbstractLinearModel. This is a base class which contains a number of methods and properties to store and compute information on a model. The AbstractLinearModel is not usable on its own as the fit_model method is not implemented. When classes such as OLSLinearModel are defined in SAILS, they inherit from AbstractLinearModel to define the helper functions before a specific fit_model method is defined. We can do the same to define a custom model fit class using an external package. We will first create a new class which inherits from AbstractLinearModel and then define a fit_model method which computes the model fit and stores the outputs in a standard form.

Here is our custom model fit class, each section is described in the comments in the code.

# Define a new class inheriting from the SAILS base model
class TSALinearModel( sails.AbstractLinearModel ):

    # Define a fit_model method using the python @classmethod decorator
    @classmethod
    def fit_model( cls, data, delay_vect):

        # Some sanity checking of the input matrix We make sure that the input
        # data is in 2d format or 3d format with a single epoch.
        # statsmodels.tsa doesn't currently support fitting multitrial data
        if data.ndim == 3:
            if data.shape[2] > 1:
                raise ValueError('This class is only implemented for single-trial data')
            # Take first trial if we have 3d data
            data = data[:,:,0]

        # Create object - classmethods act as object constructors. cls points
        # to TSALinearModel and ret is a specific instance of TSALinearModel
        # though it is currently empty.
        ret = cls()

        # The rest of this method will populate ret with information based on
        # our model fit

        # Import the model fit function
        from statsmodels.tsa.api import VAR

        # Initialise and fit model - we use a simple VAR with default options.
        # Note that we return the model and results to ret.tsa_model.  This
        # means that the statsmodels.tsa.api.VAR object will be stored and
        # returned with ret. Later we can use this to access the statsmodels
        # model and results though the SAILS object.
        ret.tsa_model = VAR(data.T) # SAILS assumes channels first, TSA samples first
        model_order = len(delay_vect) - 1 # delay_vect includes a leading zero
        ret.tsa_results = ret.tsa_model.fit(model_order)

        # The method must assign the following values for SAILS metrics to work properly
        ret.maxorder = model_order
        ret.delay_vect = np.arange(model_order)
        ret.parameters = np.concatenate((-np.eye(data.shape[0])[:,:,None],
                                          ret.tsa_results.coefs.transpose((1,2,0))), axis=2)
        ret.data_cov = sails.find_cov(data.T,data.T)
        ret.resid_cov = sails.find_cov(ret.tsa_results.resid.T,ret.tsa_results.resid.T)

        # Return fitted model within an instance of a TSALinearModel
        return ret

It is crucial that the fit_model class returns an instance of our the overall class. This instance must contain the following information. Other functions in SAILS assume that these are stored in a fitted model class with specific formats and names.

• maxorder: the model order of the fit
• delay_vect: the vector of delays used in the model fit
• parameters: the fitted autoregressive parameters of shape [num_channels x num_channels x model_order] with a leading identity
• data_cov: the covariance matrix of the fitted data
• resid_cov: the covariance matrix of the residuls of the fit

Other data can be added in as well (we store tsa_model and tsa_results in the example here) but these five must be defined within the returned class.

We can now fit a model using our new class

tsa_model = TSALinearModel.fit_model(X,np.arange(4))

Finally, we compute connectivity metrics from each model fit and plot a comparison

freq_vect = np.linspace(0,sample_rate/2)

sails_metrics = sails.FourierMvarMetrics.initialise(sails_model,sample_rate,freq_vect)
tsa_metrics = sails.FourierMvarMetrics.initialise(tsa_model,sample_rate,freq_vect)

PDC = np.concatenate( (sails_metrics.partial_directed_coherence,tsa_metrics.partial_directed_coherence),axis=3)

sails.plotting.plot_vector(PDC,freq_vect,line_labels=['SAILS','TSA'],diag=True,x_label='Frequency (Hz'))

We see that the partial directed coherence from the two models is nearly identical.

Tutorial 9 - Modal decomposition of a simulated system

One of the main features of SAILS is the ability to perform a modal decomposition on fitted MVAR models. This tutorial will outline the use of modal decomposition on some simulated data. This example is included as part of an upcoming paper from the authors of SAILS.

The system which is being simulated consists of 10 “nodes”. Each node is made up of a mixture of two components, each of which is distributed differently amongst the nodes. We start by simulating the data for a single realisation:

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

import sails

plt.style.use('ggplot')

# General configuration for the generation and analysis
sample_rate = 128
num_samples = sample_rate*50

f1 = 10
f2 = 0

siggen = sails.simulate.SailsTutorialExample()

num_sources = siggen.get_num_sources()
weight1, weight2 = siggen.get_connection_weights()

# Calculate the base signals for reference later on
sig1 = siggen.generate_basesignal(f1, .8+.05, sample_rate, num_samples)
sig2 = siggen.generate_basesignal(f2, .61, sample_rate, num_samples)

sig1 = (sig1 - sig1.mean(axis=0)) / sig1.std(axis=0)
sig2 = (sig2 - sig2.mean(axis=0)) / sig2.std(axis=0)

# Generate a network realisation
X = siggen.generate_signal(f1, f2, sample_rate, num_samples)

At this point the variable X contains the simulated data. We start by visualising this data. sig1 and sig2 contain the base signals which are used for visualistion purposes only. Finally, weight1 and weight2 contain the connection weights (the amount to which each signal is spread across the nodes).

Next we visualise the system before fitting our model. This section of the code is for plotting purposes only and can be skipped if you are happy to just examine the plot following it.

# Calculate the Welch PSD of the signals
f, pxx1 = signal.welch(sig1.T, fs=sample_rate, nperseg=128, nfft=512)
f, pxx2 = signal.welch(sig2.T, fs=sample_rate, nperseg=128, nfft=512)

# Visualise the signal and the mixing of it
x_label_pos = -2.25
x_min = -1.75
plot_seconds = 5

plot_samples = int(plot_seconds * sample_rate)

# Plot the two base signals
plt.figure(figsize=(7.5, 2.5))
plt.axes([.15, .8, .65, .15], frameon=False)
plt.plot(sig1[:plot_samples], 'r')
plt.plot(sig2[:plot_samples] + 7.5, 'b')

plt.xlim(-10, plot_samples)
yl = plt.ylim()
plt.xticks([])
plt.yticks([])

# Add labelling to the base signals
plt.axes([.05, .8, .1, .15], frameon=False)
plt.text(x_label_pos, 7.5, 'Mode 1', verticalalignment='center', color='b')
plt.text(x_label_pos, 0,   'Mode 2', verticalalignment='center', color='r')

plt.xlim(x_label_pos, 1)
plt.ylim(*yl)
plt.xticks([])
plt.yticks([])

# Add a diagram of the Welch PSD to each signal
plt.fill_between(np.linspace(0, 2, 257), np.zeros((257)), 20 * pxx1[0, :], color='r')
plt.fill_between(np.linspace(0, 2, 257), np.zeros((257)) + 7.5, 20 * pxx2[0, :]+7.5, color='b')
plt.xlim(x_min, 1)

# Plot the mixed signals in the network
plt.axes([.15, .1, .65, .7], frameon=False)
plt.plot(X[:, :plot_samples, 0].T + np.arange(num_sources) * 7.5, color='k', linewidth=.5)
plt.hlines(np.arange(num_sources)*7.5, -num_sources, plot_samples, linewidth=.2)
plt.xlim(-num_sources, plot_samples)
plt.xticks([])
plt.yticks([])
yl = plt.ylim()

# Add labelling and the weighting information
plt.axes([.05, .1, .1, .7], frameon=False)
plt.barh(np.arange(num_sources)*7.5+1, weight1, color='r', height=2)
plt.barh(np.arange(num_sources)*7.5-1, weight2, color='b', height=2)
plt.xlim(x_min, 1)
plt.ylim(*yl)
for ii in range(num_sources):
    plt.text(x_label_pos, 7.5 * ii, 'Node {0}'.format(ii+1), verticalalignment='center')
plt.xticks([])
plt.yticks([])

# Plot the mixing matrices
wmat1, wmat2 = siggen.get_connection_weights('matrix')

plt.axes [.85, .6, .13, .34])
plt.pcolormesh(wmat2, cmap='Blues')
plt.xticks(np.arange(num_sources, 0, -1) - .5, [''] * num_sources, fontsize=6)
plt.yticks(np.arange(num_sources, 0, -1) - .5, np.arange(num_sources, 0, -1), fontsize=6)
plt.axis('scaled')

plt.axes([.85, .2, .13, .34])
plt.pcolormesh(wmat1, cmap='Reds')
plt.xticks(np.arange(num_sources, 0, -1)-.5, np.arange(num_sources, 0, -1), fontsize=6)
plt.yticks(np.arange(num_sources, 0, -1)-.5, np.arange(num_sources, 0, -1), fontsize=6)
plt.axis('scaled')

Next we fit an MVAR model and examine the Modal decomposition:

model_order = 5
delay_vect = np.arange(model_order + 1)
freq_vect = np.linspace(0, sample_rate/2, 512)

# Fit our MVAR model and perform the modal decomposition
m = sails.VieiraMorfLinearModel.fit_model(X, delay_vect)

modes = sails.MvarModalDecomposition.initialise(m, sample_rate)

# Now extract the metrics for the modal model
M = sails.ModalMvarMetrics.initialise(m, sample_rate,
                                      freq_vect, sum_modes=False)

Each mode in the decomposition consists of either a single pole (where it is a DC component) or a pair of poles (where it is an oscillator). We can extract the indicies of the pole pairs using the `mode_indices` property. Given the nature of our simulation, we would expect to find one DC pole and one pole at 10Hz with a stronger effect than all other poles. For the purposes of this tutorial, we are going to set an arbitrary threshold on the modes; in reality a permutation testing scheme can be used to determine this threshold (which will be the subject of a future tutorial).

# For real data and simulations this can be estimated using
# a permutation scheme
pole_threshold = 2.3

# Extract the pairs of mode indices
pole_idx = modes.mode_indices

# Find which modes pass threshold
surviving_modes = []

for idx, poles in enumerate(pole_idx):
   # The DTs of a pair of poles will be identical
   # so we can just check the first one
   if modes.dampening_time[poles[0]] > pole_threshold:
      surviving_modes.append(idx)

# Use root plot to plot all modes and then replot
# the surviving modes on top of them
ax = sails.plotting.root_plot(modes.evals)

low_mode = None
high_mode = None

for mode in surviving_modes:
   # Pick the colour based on the peak frequency
   if modes.peak_frequency[pole_idx[mode][0]] < 0.001:
       color = 'b'
       low_mode = mode
   else:
       color = 'r'
       high_mode = mode

   for poleidx in pole_idx[mode]:
      ax.plot(modes.evals[poleidx].real, modes.evals[poleidx].imag,
              marker='+', color=color)

From this, we can see that the modal decomposition has clearly extracted a set of poles of importance. We can now visualise the connectivity patterns for these modes as well as their spectra:

# The frequency location simply acts as a scaling factor at this
# point, but we pick the closest bin in freq_vect to the
# peak frequency of the mode

low_freq_idx = np.argmin(np.abs(freq_vect - modes.peak_frequency[pole_idx[low_mode][0]]))
high_freq_idx = np.argmin(np.abs(freq_vect - modes.peak_frequency[pole_idx[high_mode][0]]))

# We can now plot the two graph
plt.figure()

low_psd = np.sum(M.PSD[:, :, :, pole_idx[low_mode]], axis=3)
high_psd = np.sum(M.PSD[:, :, :, pole_idx[high_mode]], axis=3)

# Plot the connectivity patterns as well as the spectra
plt.subplot(2, 2, 1)
plt.pcolormesh(low_psd[:, :, low_freq_idx], cmap='Blues')
plt.xticks(np.arange(num_sources, 0, -1)-.5, np.arange(num_sources, 0, -1), fontsize=6)
plt.yticks(np.arange(num_sources, 0, -1)-.5, np.arange(num_sources, 0, -1), fontsize=6)

plt.subplot(2, 2, 2)
for k in range(num_sources):
   plt.plot(freq_vect, low_psd[k, k, :])

plt.subplot(2, 2, 3)
plt.pcolormesh(high_psd[:, :, high_freq_idx], cmap='Reds')
plt.xticks(np.arange(num_sources, 0, -1)-.5, np.arange(num_sources, 0, -1), fontsize=6)
plt.yticks(np.arange(num_sources, 0, -1)-.5, np.arange(num_sources, 0, -1), fontsize=6)

plt.subplot(2, 2, 4)
for k in range(num_sources):
   plt.plot(freq_vect, high_psd[k, k, :])

plt.xlabel('Frequency (Hz)')

In this tutorial we have demonstrated how to use the modal decomposition functionality in SAILS to decompose the parameters of an MVAR model and then to extract spectral and connectivity patterns from just the modes of significance. This approach can be extended across participants and, with the additional use of PCA, used to simultaneously estimate spatial and spectral patterns in data. More details of this approach will be available in an upcoming paper from the authors of SAILS.

Tutorial 10 - Dynamic connectivity during a task

Here, we will look at using MVAR modelling to describe changes in connectivity within a functional network as participants perform a simple button press task. This is similar to the sliding window modelling in tutorial 6

We will analyse MEG source time-courses from four regions of the AAL atlas (Precentral gyrus and supplemental motor area from left and right hemispheres) during a self-paced finger tap task from 10 participants. Each trial lasts 20 seconds with 10 seconds of finger tapping at the start and 10 seconds post movement time. Finger tapping was performed with the right hand. The MEG data were recorded on a 4D NeuroImaging WHS-3600 system and the source time-courses were generated from the data using an LCMV beamformer on data which had been band-pass filtered between 1 and 80Hz.

First, lets import sails and load in the example data. If you haven’t already done so, please download the example data repository from https://vcs.ynic.york.ac.uk/analysis/sails-example-data

We start by importing the modules we will require:

import os

import h5py
import numpy as np
import matplotlib.pyplot as plt

import sails

SAILS will automatically detect the example data if downloaded into your home directory. If you’ve used a different location, you can specify this in an environment variable named SAILS_EXAMPLE_DATA.

# Specify environment variable with example data location
# This is only necessary if you have not checked out the
# example data into $HOME/sails-example-data
os.environ['SAILS_EXAMPLE_DATA'] = '/path/to/sails-example-data'

# Locate and validate example data directory
example_path = sails.find_example_path()

# Load data using h5py
motor_data = h5py.File(os.path.join(sails.find_example_path(), 'fingertap_group_data.hdf5'))

The motor data is stored in hdf5 format and contains the data sample rate and 10 data arrays with the data for each participant. These can be accessed using keys similar to a dictionary. Here, we print the keys from motor_data and extract the sample_rate. Note that the motor_data['sample_rate'] returns a h5py object which we can further index to extract the sample rate using motor_data['sample_rate'][0]

# Print contents of motor_data
print(list(motor_data.keys()))

# Extract sample_rate
sample_rate = motor_data['sample_rate'][...]
print('Data sample rate is {0}Hz'.format(sample_rate))

# Define node labels
labels = ['L Precentral', 'R Precentral', 'L SuppMotorArea', 'R SuppMotorArea']

The fingertap data itself is in a 3d array of size [nchannels x nsamples x ntrials]. Every participant has 4 channels and 3391 samples in each trial but slightly different numbers of trials - around 20-30 each.

# Print shape of data array from the first participant
print(motor_data['subj0'][...].shape)

Before fitting our model we specify a time vector with the time in seconds of each of our samples.

# Specify a time vector
num_samples = motor_data['subj0'].shape[1]
time_vect = np.linspace(0, num_samples/sample_rate, num_samples)

Now we will fit our models. We first define the vector of delays to fit the MVAR model on and a set of frequency values to estimate connectivity across. We will compute three things for each participant: m is the LinearModel containing the autoregressive parameters, d is a set of model diagnostics for each mode and f is a MvarMetrics instance which we can use to compute power and connectivity values.

We compute m, d and f for each participant in turn and store them in a list. Please see tutorial 6 for more details on sliding_window_fit and its options.

# Define model delays, time vector and frequency vector
delay_vect = np.arange(15)
freq_vect = np.linspace(0, 48, 48)

# Initialise output lists
M = []
D = []
F = []

# Main loop over 10 subjects
for ii in range(10):
    print('Processing subj {0}'.format(ii))

    # Get subject data
    x = motor_data['subj{}'.format(ii)][...]

    # Fit sliding window model
    sliding_window_length = int(sample_rate) # 1 second long windows
    sliding_window_step = int(sample_rate / 8) # 125ms steps between windows
    m, d = sails.sliding_window_fit(sails.VieiraMorfLinearModel, x, delay_vect,
                                    sliding_window_length, sliding_window_step)

    # Compute Fourier MVAR metrics from sliding window model
    f = sails.FourierMvarMetrics.initialise(m, sample_rate, freq_vect)

    # Append results into list
    M.append(m) # MVAR Model
    D.append(d) # Model Diagnostics
    F.append(f) # Fourier Metrics

# Get time vector for centre of sliding windows (in seconds)
model_time_vect = time_vect[m.time_vect.astype(int)]

We can extract information across participants using list comprehensions. Here, we extract the power spectral density from each participant and concatenate them into a single array for visualisation.

# Create a list of PSD arrays with a singleton dummy dimension on the end
# and concatenate into a single array
PSD = np.concatenate([x.PSD[..., np.newaxis] for x in F], axis=4)

# PSD is now [nnodes x nnodes x nfrequencies x ntimes x nparticipants]
print(PSD.shape)

Next we visualise the time-frequency power spectral density for each of the four nodes. We perform a simple baseline correction by subtracting the average of the last 2 seconds of data from the whole trial. The resulting PSD shows the power relative to this pre-movement period. We annotate the plots with two dotted lines, one at 10 seconds to show the end of the finger-tapping and one at 18 seconds showing the start of the baseline period.

# Count the number of nodes and subjects
num_nodes = PSD.shape[0]

# Number of windows over which to calculate baseline estimate
baseline_windows = 11

# Create a new figure
plt.figure(figsize=(6, 12))

# Main plotting loop
for ii in range(num_nodes):
    # Average PSD across  participants
    psd = PSD[ii, ii, :, :, :].mean(axis=2)

    # Apply a simple baseline correction
    psd = psd - psd[:, -baseline_windows:, np.newaxis].mean(axis=1)

    # Find the max value for the colour scale
    mx = np.abs(psd).max()

    # Make new subplot and plot baseline corrected PSD
    plt.subplot(num_nodes, 1, ii + 1)
    plt.pcolormesh(model_time_vect, freq_vect, psd, cmap='RdBu_r', vmin=-mx, vmax=mx)

    # Annotate subplot
    cb = plt.colorbar()
    cb.set_label('PSD')

    # Place lines showing the period of finger tapping
    plt.vlines([10, 18], freq_vect[0], freq_vect[-1], linestyles='dashed')

    # Annotate windows
    plt.text(5, 40, 'Tapping', horizontalalignment='center')
    plt.text(14, 40, 'Rebound', horizontalalignment='center')
    plt.text(18.75, 40, 'BL', horizontalalignment='center')

    # Tidy up x-axis labelling
    if ii == (num_nodes - 1):
        plt.xlabel('Time (seconds)')
    else:
        plt.gca().set_xticklabels([])

    # Y axis labelling and title
    plt.ylabel('Frequency (Hz)')
    plt.title(labels[ii])

plt.show()

Note that the Left Precentral gyrus has a strong increase in beta power after movement has stopped. The left and right Supplemental Motor Areas have a weaker rebound.

It is always good idea to inspect the model diagnostic values for an MVAR analysis. We now extract the stability index, r-squared and residual covariances for each participant using list comprehensions to extract data from D.

We use the np._r operator as a quick way to concatenate our lists into numpy arrays.

# Get stability index
SI         = np.r_[[d.SI for d in D]]

# Get R-square variance explained
R_square   = np.r_[[d.R_square.mean(axis=1) for d in D]]

# Get the matrix norm of the residual covariance matrices - this is a
# convenient summary of the sum-squared values in the residual covariance
# matrices.
resid_norm = np.r_[[np.linalg.norm(d.resid_cov, axis=(0, 1)) for d in D]]

A quick visualisation of these diagnostics shows that our models are stable for all participants and all time windows (SI < 1). The models explain between 15 and 40% of variance and have relatively stable residual covariances across the whole window.

plt.figure()

plt.subplot(3, 1, 1)
plt.plot(model_time_vect, SI.T, 'grey')
plt.plot(model_time_vect, SI.mean(axis=0), 'k', linewidth=2)
plt.ylabel('Stability Index')
plt.gca().set_xticklabels([])

plt.subplot(3, 1, 2)
plt.plot(model_time_vect, R_square.T, 'grey')
plt.plot(model_time_vect, R_square.mean(axis=0), 'k', linewidth=2)
plt.ylabel('R-square')
plt.gca().set_xticklabels([])

plt.subplot(3, 1, 3)
plt.plot(model_time_vect, resid_norm.T, 'grey')
plt.plot(model_time_vect, resid_norm.mean(axis=0), 'k', linewidth=2)
plt.ylabel('Norm of\nresidual covariance')

plt.show()

Now we trust that our models are capturing reasonable task dynamics within each brain region and have good diagnostics we can look at the connectivity.

We first look at the cross-spectral densities across the network. These are the off diagonal elements of the PSD metric. We first extract the PSD using the list comprehension method and concatenate them into a single array. After that, we plot the average cross spectral density for between all nodes using sails.plotting.plot_matrix.

# Create a list of PSD arrays with a singleton dummy dimension on the end
# and convert into an array
PSD = np.concatenate([f.PSD[..., np.newaxis] for f in F], axis=4)

# Visualise
fig = plt.figure(figsize=(12, 8))
sails.plotting.plot_matrix(PSD.mean(axis=4), model_time_vect, freq_vect,
                           title='Cross Spectral Density',
                           labels=labels, F=fig,
                           vlines=[10], cmap='hot_r', diag=False,
                           x_label='Time (secs)', y_label='Frequency (Hz)')
fig.show()

The cross spectral densities show a similar post-movement beta rebound pattern to the within-node power spectral densities. Now we can also see that there is shared spectral information in the left-precentral gyrus <-> left-supplemental motor area and left-supplemental motor area <-> right-supplemental motor area connections. There appears to be strong cross-spectral densities below 10Hz between all nodes.

The Magnitude-Squared Coherence might be a better representation of these connections. It expresses the cross-spectral density between two nodes as a ratio of the power within each node.

# Extract the magnitude squared coherence using the list comprehension method
# and convert into a numpy array
MSC = np.concatenate([f.magnitude_squared_coherence[..., np.newaxis] for f in F], axis=4)

# Visualise
fig = plt.figure(figsize=(12, 8))
sails.plotting.plot_matrix(MSC.mean(axis=4), model_time_vect, freq_vect,
                           title='Magnitude Squared Coherence',
                           labels=labels, F=fig,
                           vlines=[10], cmap='hot_r', diag=False,
                           x_label='Time (secs)', y_label='Frequency (Hz)')
plt.show()

The normalisation emphasises the coherence within the beta rebound and strongly reduces the apparent shared power below 10Hz. This suggests that the beta cross spectral density is relatively large when compared to the power in each node at that frequency, but the <10Hz cross spectra are very low power compared to the within node power.

Next, we can explore whether this beta connectivity is symmetrical i.e. whether both nodes are equally influential on each other or if one node in the pair might be ‘driving’ the other. We use the Directed Transfer Function to estimate this and visualise in the same way.

# Extract the directed transfer function using the list comprehension method
# and convert into a numpy array
DTF = np.concatenate([f.directed_transfer_function[..., np.newaxis] for f in F], axis=4)

# Visualise
fig = plt.figure(figsize=(12, 8))
sails.plotting.plot_matrix(DTF.mean(axis=4), model_time_vect, freq_vect,
                           title='Directed Transfer Function',
                           labels=labels, F=fig,
                           vlines=[10], cmap='hot_r', diag=False,
                           x_label='Time (secs)', y_label='Frequency (Hz)')
plt.show()

The DTF is an asymmetrical measure, so the upper and lower triangles of the DTF plot are not symmetrical. We see similar connections in the beta band again, but the DTF additionally suggests that Left Precentral Gyrus which is driving Left Supplemental Motor Area, though there is some influence in the reciprocal direction. Similarly Left Supplemental Motor Area appears to be influencing Right Supplemental Motor Area.

Finally, we can emphasise the change in connectivity relative to baseline by performing a simple baseline correction on the DTF values. Here, we subtract the average DTF from the last two seconds of the epoch from each time-point. Positive values then indicate a movement-evoked increase in connectivity in a connection and negative values a movement-evoked decrease.

# Number of windows over which to calculate baseline estimate
baseline_windows = 11

# Apply a simple baseline correction
bcDTF = DTF.mean(axis=4)
bcDTF = bcDTF - bcDTF[:, :, :, -baseline_windows:, np.newaxis].mean(axis=3)

# Plot baseline corrected DTF
fig = plt.figure(figsize=(12, 8))
sails.plotting.plot_matrix(bcDTF, model_time_vect, freq_vect,
                           title='baseline-corrected Directed Transfer Function',
                           labels=labels, F=fig,
                           vlines=[10, 18], cmap='RdBu_r', diag=False,
                           x_label='Time (secs)', y_label='Frequency (Hz)')
plt.show()

This baseline correction makes the change in directed functional connectivity during the post-movement beta rebound much clearer. It also reveals the fact that the relationship between the two supplementary motor areas appears to be driven by the left SMA. Given that this is a right-hand movement task, this could potentially be interpreted as a form of inhibitory signal from the left to the right hemisphere. Further data and analysis would be necessary to fully establish the nature of such a signal.

Gallery

Simulation Gallery

This page contains a set of figures re-created around existing papers which used simulated or analytic systems. These graphs are provided for comparison with the original papers.

The code for generating these figures can be found in sim_gallery.py.

Baccala and Sameshima 2001

We re-create figures 1 to 5 of this paper. Each figure compares PDC and DTF connectivity estimators.

Pascual-Marqui et al 2014

We re-create figure 3 of this paper using IEC and PDC. PDC and DTF connectivity estimators.

API

The sails API can be divided up into a number of categories, each of which has its own page:

Data Simulation

class sails.simulate.AbstractSigGen(extraargs=None)

Abstract base class for signal generation

generate_signal(meg_reader)

Generate a custom signal for embedding in a dataset

Parameters:meg_reader – MEG data for which to generate signal Returns:numpy array of shape (samples). Note that samples must be equal to the value of meg_reader.num_samples Return type:type

get_num_sources()

Helper function for returning number of sources in simulation.

This is typically overloaded when defining a signal generator class.

Returns:type – Number of sources in simulated signal. If not override, returns 1 Return type:int

parse_arguments(extraargs)

Helper function for processing extra-arguments passed to a signal-generator.

Currently only looks for snr_dB for noise levels but can be overloaded by a specific class to processes more arguments.

Parameters:extraargs (dict) – keyword arguments passed to signal generator.

class sails.simulate.AbstractLagSigGen(extraargs=None)

Generic class for making lagged signals

generate_basesignal(nsamples)

Function defining the starting signal which will be lagged throughout the network

This is an abstracted function and should be overloaded when creating a new class.

Parameters:nsamples (int) – The number of samples to generate

generate_noise(nsamples)

Function defining the additive noise added at each lag

Parameters:nsamples (int) – The number of samples to generate Returns:A vector of white noise Return type:ndarray

generate_signal(sample_rate=None, num_samples=None, noise_dB=None)

Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class

Returns a timeseries of dimension [samples x signals]

Parameters:

• sample_rate (float) – The sample rate to generate the signal at (Default value = None)
• num_samples – The number of samples to generate (Default value = None)
• noise_dB (float) – Noise level of the additive noise (Default value = None)

Returns:

A simulated time-series [nchannels x nsamples x nrealisations]

Return type:

ndarray

get_lags()

Method returning the lagging of the signal through the system. Should

Returns:dimension refers to the source node and the second dimension the target node. Each non-zero value indicates that the signal in the first dimension will be added into the signal in the second dimension at the specified weight. Return type:type

get_num_sources()

parse_arguments(extraargs)

Helper function for processing extra-arguments passed to a signal-generator.

Currently only looks for snr_dB for noise levels but can be overloaded by a specific class to processes more arguments.

Parameters:extraargs (dict) – keyword arguments passed to signal generator.

class sails.simulate.AbstractMVARSigGen(extraargs=None)

Generic class for making MVAR signals

check_stability()

Checks the stability of the parameter matrix by the magnitude of the largest eigenvalue of the first lag in the parameter matrix (see [Lutkephol2006] chapter 1)

Returns:The stability index of the simulated signal Return type:float

Notes

The stability index of a simulation should be strictly less than 1. If this exceeds 1 then the simulated signal is likely to tend to infinity or zero.

generate_model()

Returns a LinearModel containing containing the true model parameters.

generate_parameters()

Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

This is an abstract method and should be overloaded when creating a specific class.

generate_signal(sample_rate=None, num_samples=None, noise_dB=None, num_realisations=1)

Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class

Returns a timeseries of dimension [nchannels x nsamples x nrealisations]

Parameters:

• sample_rate (float) – The sample rate of the simulated data (Default value = None)
• num_samples (int) – The number of samples to generate (Default value = None)
• noise_dB (float) – The level of additive noise to add to the simulationo (Default value = None)
• num_realisations (int) – Number of realisations to generate (Default value = 1)

Returns:

The simulated dataset of size [nchannels x nsamples x nrealisations]

Return type:

ndarray

get_num_sources()

parse_arguments(extraargs)

Helper function for processing extra-arguments passed to a signal-generator.

Currently only looks for snr_dB for noise levels but can be overloaded by a specific class to processes more arguments.

Parameters:extraargs (dict) – keyword arguments passed to signal generator.

class sails.simulate.Baccala2001_fig1(extraargs=None)

This class implements [Baccala2001] figure 1.

generate_parameters()

Create the autoregressive parameters for [Baccala2001] figure 1.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Baccala2001_fig2(extraargs=None)

This class implements [Baccala2001] figure 2.

generate_parameters()

Create the autoregressive parameters for [Baccala2001] figure 2.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Baccala2001_fig3(extraargs=None)

This class implements [Baccala2001] figure 3.

generate_parameters()

Create the autoregressive parameters for [Baccala2001] figure 3.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Baccala2001_fig4(extraargs=None)

This class implements [Baccala2001] figure 4.

generate_parameters()

Create the autoregressive parameters for [Baccala2001] figure 4.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Baccala2001_fig5(extraargs=None)

This class implements [Baccala2001] figure 5.

generate_parameters()

Create the autoregressive parameters for [Baccala2001] figure 5.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Fasoula2013_eqn26(extraargs=None)

This class implements [Fasoula2013] equation 26.

generate_parameters()

Create the autoregressive parameters for [Fasoula2013] eqn 26.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Schelter2006_fig1(extraargs=None)

This class implements [Schelter2006] figure 1.

generate_parameters()

Create the autoregressive parameters for [Schelter2006] fig 1.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.PascualMarqui2014_fig3(extraargs=None)

This class implements [Pascual-Marqui2014] figure 3.

generate_parameters()

Create the autoregressive parameters for [Pascual-Marqui2014] fig 3.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.Korzeniewska2003Lag(extraargs=None)

Generate realisations of the time-lagged network defined in [Korzeniewska2003] figure 1.

generate_basesignal(nsamples)

Generate oscillatory signal to be propagated through the network. This creates a single resonance at one-half of the Nyquist frequency via direct pole placement.

Note: the original base-signal in the paper is a realisation from an AR model. I couldn’t find the parameters for this mode, so use pole placement here. The spectral profile will be different to the paper, but the connectivity patterns are the same.

Parameters:nsamples (int) – The number of samples to generate Returns:Vector containing generated base signal. Return type:ndarray

get_lags()

Return the lagged influences between time-series in the simulated defined in [Korzeniewska2003] figure 1.

get_num_sources()

Return the number of sources in the simulation

class sails.simulate.SailsTutorialExample(extraargs=None)

generate_basesignal(f, r, sample_rate, num_samples)

Generate a simple signal by pole placement

Parameters:

• f (float) – The peak frequency of the resonance to simulate
• r (float ( 0 < r < 1 )) – The pole magnitude of the resonance to simulate
• sample_rate (float) – The sample rate to generate data at
• num_samples (int) – The number of data samples to generate

Returns:

A vector containing a simulated resonance

Return type:

ndarray

Notes

The parameter r controls both the sharpness and magnitude of the generated resonance. Values close to 1 will generate exponentially sharper and larger peaks.

generate_signal(f1, f2, sample_rate, num_samples, noise_dB=None, num_realisations=1, magnitude_noise=0.02)

Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class

Returns a timeseries of dimension [samples x signals x realisations]

Parameters:

• f1 (float) – The frequency of the first signal
• f2 (float) – The frequency of the second signal
• sample_rate (float) – The sample rate of the simulated signal
• num_samples (int) – The number of samples to simulate
• num_realisations (int) – The number of realisations to simulate (Default value = 1)
• magnitude_noise (float) – The variance of additive white noise (Default value = .02)

Returns:

The simulated signal of size [nchannels x num_samples x num_realisations]

Return type:

ndarray

get_connection_weights(form='vector')

Define the connectivity within the simulated system.

Parameters:form ({'vector' || 'matrix'}) – Which form of to return the connections is (Default value = ‘vector’) Returns:

• weight1 – vector or matrix defining connectivity for resonance 1
• weight2 – vector or matrix defining connectivity for resonance 2

get_num_sources()

Return the number of sources in the simulation

Preprocessing Functions

class sails.utils.PCA(data=None, npcs=None)

Class for handling PCA-based reduction before fitting a model

project_score(scores)

Compute a projection of a set of scores back to original data space

Model Fitting

Model fitting is normally performed by either the sails.modelfit.VieiraMorfLinearModel or sails.modelfit.OLSLinearModel class. If you wish to fit a model using another method, an example of this can be found in :ref:tutorial8.

class sails.modelfit.VieiraMorfLinearModel

A class implementing the Vieira-Morf linear model fit

classmethod fit_model(data, delay_vect)

Estimates the multichannel autoregressive spectral estimators using a Vieira-Morf algorithm. Equations are referenced to [Marple1987], appendix 15.B.

This is the multitrial versions of the algorithm, using the AMVAR method outlined in [Ding2000].

Parameters:

• data (numpy.ndarray) – array of shape [nsignals, nsamples, ntrials]
• delay_vect (numpy.ndarray) – Vector containing evenly spaced delays to be assessed in the model. Delays are represented in samples. Must start with zero.

Returns:

A populated object containing the fitted forward and backwards coefficients and several other useful variables and methods.

Return type:

sails.VieiraMorfLinearModel

get_residuals(data, forward_parameters=True, mode='valid')

Returns the prediction error from a fitted model. This is a wrapper function for get_residuals()

Parameters:

• data (ndarray) – Data to compute the residuals from, of size [nchannels x nsamples x nrealisations]
• forward_parameters (bool) – If True, use forward parameters, otherwise use backward parameters (Default value = True)

Returns:

Residual data

Return type:

ndarray

class sails.modelfit.OLSLinearModel

A class implementing ordinary least squares linear model fit

classmethod fit_model(data, delay_vect, estimator=None)

This is a class method which fits a linear model and returns a populated OLSLinearModel instance containing the fitted model.

Parameters:

• data (ndarray) – The data to be used to compute the model fit
• delay_vect (ndarray) – A vector of lag indices defining the lags to fit
• estimator (None or sklearn class) – If None, fit using standard OLS normal equations. If set to an appropriate sklearn class, use that estimator.

Returns:

A populated object containing the fitted coefficients and several other useful variables and methods.

Return type:

sails.OLSLinearModel

Model fitting helper functions

In addition to basic model fitting routines, some helper functions exist; for instance to help with fitting a number of models to a data series following a “sliding window” pattern or by applying PCA before the fit.

sails.modelfit.sliding_window_fit(model_class, data, delay_vect, win_len_samp, win_step_samp, compute_diagnostics=True)

A helper function for fitting many MVAR models within sliding windows across a dataset.

Parameters:

• model_class (LinearModel) – The SAILS linear model typee to fit
• data (ndarray) – Data to fit the model on, of size [nchannels x nsamples x nrealisations]
• delay_vect (ndarray) – A vector of lags specifying which lags to fit
• win_len_samp (int) – The window length in samples
• win_step_samp (int) – The step size between windows in samples
• compute_diagnostics (boolean) – Flag indicating whether to compute model diagnostics at each window (Default value = True)

Returns:

An instance of linear model passed into the function containing MVAR parameters for all windows. The parameters are stored in model.parameters which is of size: [nchannels x nchannels x model order x nwindows]

Return type:

LinearModel

sails.modelfit.pca_reduced_fit(X, delay_vect, ndim, linear_model=<class 'sails.modelfit.VieiraMorfLinearModel'>)

Helper for computing an MVAR on dimensionality reduced data (using PCA).

Returns a model fitted to reduced data, the reduced model projected back to original data dimensions and the pca object used for reduction

Parameters:

• X (ndarray) – The data to compute the reduced model fit on
• delay_vect (ndarray) – A vector of lags specifying which lags to fit
• ndim (int) – Number of components to reduce data to.
• linear_model (class) – Subclass of AbstractLinearModel to use for fit. Defaults to VieiraMorfLinearModel. (Default value = True)

Returns:

• red_model (AbstractLinearModel) – An instance of the linear model passed into the function containing MVAR parameters for all windows for the reduced model.
• proj_model (AbstractLinearModel) – An instance of the linear model passed into the function containing MVAR parameters for all windows for the projected model.
• pc (sails.utils.PCA) – PCA class with information about the PCA projection included

Model Diagnostics

Several classes and/or routines are available to help with examining and diagnosing model fits.

class sails.diags.ModelDiagnostics

Class to store and display LinearModel diagnostic information. Several diagnostic criteria are computed including:

R_square (R_square)

the percentage of variance explained in each channel

Stability Index (SI)

Indicator of the stability of the MVAR parameters (SI<1 indicates a stable model)

Stability Ratio (SR)

A stronger test of stability computed from the ratio of the largest eigenvalue of A to all others.

Durbin-Watson (DW)

A test of autocorrelation in residuals of the model fit. Values should be close to 2 indicating no autocorrelation, values close to 0 indicate a positive autocorrelation and 4 and negative autocorrelation.

Log Likelihood (LL)

The log-likelihood of the model.

Akaike’s Information Criterion (AIC)

An indication of the model ‘quality’, lower values indicate a more accurate, less complex model.

Bayesian Information Criterion (BIC)

An indication of the model ‘quality’, lower values indicate a more accurate, less complex model.

Percent Consistency (PC)

Indicates how well a model captures the auto and cross correlation in a time-series. Only computed if compute_pc is passed to the relevant function.

classmethod combine_diag_list(diags)

Helper function for combining diagnostics from a list of ModelDiagnostics instances for easy comparison and visualisation.

Parameters:diags (list of ModelDiagnostics instances) – The ModelDiagnostics to concatenate Returns: Return type:sails ModelDiagnostics instance

classmethod compute(model, data, compute_pc=False)

Classmethod for computing a set of model diagnostics from a fitted model applied to a time-series dataset.

Parameters:

• model (sails LinearModel class) – A fitted linear model
• data (ndarray) – A 3d time-series of size [nchannels x nsamples x ntrials]
• compute_pc (bool) – Flag indicating whether to compute the percent consistency, this can be time-consuming for large datasets (Default=False).

Returns:

Return type:

sails ModelDiagnostics instance

summary(all_models=True)

Print the ModelDiagnostics in a pre-formatted table

class sails.diags.DelayDiagnostics

Class which computes the mutual information as a function of lag from zero lag to the first zero crossing of the autocorrelation function.

classmethod delay_search(data, maxdelay, step, sample_rate, constant_window=True)

Compute MI as a function of lag from zero lag to the first zero crossing of the autocorrelation function

Parameters:

• data (numpy.ndarray) – array of [signals, samples, trials]
• maxdelay (int) – maximum delay to consider
• step (int) – step to increment the delay by
• sample_rate (float) – sample rate of data
• constant_window (bool) – Flag indicating that the same number of datapoints should be included at each delay. Default is True

Returns:

Populated object containing diagnostics for each value in delay

Return type:

DelayDiagnostics

sails.modelfit.get_residuals(data, parameters, delay_vect, backwards=False, mode='valid')

This is a helper function for computing the residuals of a dataset after the MVAR predictions have been removed.

Parameters:

• data (ndarray) – Data to compute the residuals from, of size [nchannels x nsampes x nrealisations]
• parameters (ndarray) – MVAR parameter matrix, of size [nchannels x nchannels x model order]
• delay_vect (ndarray) – Vector of lag indices corresponding to the third dimension of the parameter matrix
• backwards (bool) – Flag indicating whether the forwards or backwards parameters havebeeen passed (Default value = False)
• mode ({'valid','full_nan','full'}) – Options for excluding or replacing residuals which do not have a full model prediction ie the third sample of an order 5 model. ‘valid’ removes samples without a full model prediction, ‘full_nan’ returns resids of the same size as the data with nans replacing excluded samples and ‘full’ returns resids keeping non-full model samples in place.

Returns:

Residual data, of size [nchannels x nsamples x nrealisations]

Return type:

ndarray

Model Decomposition

class sails.modal.MvarModalDecomposition

Object which calculates the modal decomposition of a fitted model derived from the AbstractLinearModel class. This is typically created using the initialise() classmethod which computes the decomposition and returns a populated object instance.

compute_modal_parameters()

Compute the order 1 or 2 autoregressive coefficients associated with each pole or conjugate pair.

Warning

This function is a work in progress and has not been fully validated.

Returns:Array containing the time-domain AR parameters for each mode Return type:ndarray

compute_residues()

Compute the residue matrix from the eigenvectors associated with each pole.

Returns:Array containing residue matrices for each pole. Return type:ndarray

excitation(resid_cov)

Compute the mode excitation as defined in [Neumaier2001]. This is a metric to quantify the dynamical importance of each mode in the decomposition.

Warning

This function is a work in progress and has not been fully validated.

Parameters:resid_cov (ndarray) – Residual covariance matrix for modelled system Returns:Vector containing modal excitation values Return type:ndarray

get_mode_inds(fmin=0, fmax=None, mag_thresh=0, resid_thresh=0, index_mode='inclusive')

Helper function to identify mode indices based on given criteria.

Parameters:

• fmin (float) – Minimum frequency of modes to include (Default value = 0)
• fmax (float) – Maximum frequency of modes to include (Default value = None)
• mag_thresh (float ( 0 > mag_thresh > 1 )) – Minimum pole magnitude to include (Default value = 0)
• resid_thresh (float) – Minimum value of mode residue norm to include (Default value = 0)
• index_mode ({'inclusive','exclusive'}) – Flag indicating whether mode selection limits should be inclusive or exclusive (Default value = ‘inclusive’)

Returns:

Array of integers indexing the included poles

Return type:

ndarray

classmethod initialise(model, sample_rate=None, normalise=False)

Compute the modal pole-residue decomposition of the A matrix from a fitted MVAR model.

Currently only implemented for a single realisation (A.ndim == 3)

Parameters:

• model (LinearModel) – Fitted model object. Must be derived from AbstractLinearModel.
• sample_rate (float) – sample rate of system in Hz (Default value = None)
• normalise (bool) – Whether to adjust the phase of the eigenvectors (Default value = False)

Returns:

Modal decomposition object

Return type:

type

modal_freqz()

Compute filter characteristics for each pole

modal_timecourse(data)

Compute the modal time-course from a dataset and a fitted model. This is the transformation of a set of [nsignals x time] data observations into modal coordinates of shape [nmodes x time].

Warning

This function is a work in progress and has not been fully validated.

Parameters:data (ndarray) – Input data to convert into modal co-ordinates Returns:Data transformed into modal time-series Return type:ndarray

modal_transfer_function(sample_rate, modes=None)

Compute a transfer function matrix based on the excitation period for each mode

Parameters:

• sample_rate (float) – sample rate of system in Hz
• modes (list of int) – List of mode indices to evaluate over (optional) (Default value = None)

Returns:

Transfer function computed for each mode

Return type:

ndarray

mode_indices

Returns a list of tuples of mode indices where each tuple contains the indices of either a pole-pair or an unpaired pole

Returns:list of tuples containing indices into the modes Return type:List

per_mode_transfer_function(sample_rate, freq_vect)

Extracts the transfer function for each mode by summing the transfer function across pole-pairs where necessary

The transfer function is computed for each pole using transfer_function() before summing individual modes (real valued poles or complex-conjugate pairs).

Parameters:

• sample_rate – sample rate of system in Hz
• freq_vect – Vector of frequencies at which to evaluate function

Returns:

Array containing the transfer function for each individual mode

Return type:

ndarray

period_hz

This is an old and deprecated way of accessing peak_frequency

pole_plot(ax=None, normscale=50, plottype='unitcircle')

Plot the poles of the current system.

Parameters:

• ax (matplotlib axes handle) – Optional axis on which to plot (Default value = None)
• normscale (float) – Scaling factor for points on plot (Default value = 50)
• plottype ({'unitcircle','eigenspectrum'}) – Flag indicating whether to plot results on a unit-circle or an eigenspectrum (Default value = ‘unitcircle’)

Returns:

Reference to axes on which plot was drawn

Return type:

matplotlib axes handle

residue_norm

The matrix-norm of the residue matrix of each mode

transfer_function(sample_rate, freq_vect, modes=None, sum_modes=True)

Compute the transfer function in pole-residue form, splitting the system into modes with separate transfer functions. The full system transfer function is then a linear sum of each modal transfer function.

When sum_modes is False, this function returns the transfer function for each individual pole (ie will return a transfer function for each pole in a conjugate pair). Please use per_mode_transfer_function to compute the transfer function for individual modes (ie one transfer function for each real pole or complex-conjugate pair)

Parameters:

• sample_rate (float) – sample rate of system in Hz
• freq_vect (ndarray) – Vector of frequencies at which to evaluate function
• modes (list of ints) – List of mode indices to evaluate over (optional) (Default value = None)
• sum_modes (bool) – Boolean indicating whether to sum modes or return all (Default value = True)

Returns:

Array containing the transfer function for each individual pole

Return type:

ndarray

Model metrics

class sails.mvar_metrics.FourierMvarMetrics

This class for computing various connectivity metrics based on the Fourier transform of the MVAR co-efficients.

This is typically called using the initialise() classmethod to compute the frequency transform and return a class instance with methods for computing different connectivity estimators

classmethod initialise(model, sample_rate, freq_vect, nmodes=None)

Compute some basic values from the Fourier transform of the A matrix. This is a class method which creates and returns a FourierMvarMetrics instance with methods computing various connectivity metrics.

This class computes the transfer function H the Fourier transform.

Parameters:

• model (sails LinearModel) – SAILS class instance containing a fitted linear model.
• sample_rate (float) – Sample rate of the dataset
• freq_vect (ndarray) – Vector of specifying which frequencies to compute connectivity at.

Returns:

Return type:

sails.FourierMvarMetrics instance

class sails.mvar_metrics.ModalMvarMetrics

This class for computing various connectivity metrics based on a Modal decomposition of a fitted MVAR model

This is typically called using the initialise() or initialise_from_modes() classmethods to compute the decomposition and return a class instance containing them.

classmethod initialise(model, sample_rate, freq_vect, sum_modes=True)

Compute some basic values from the pole-residue decomposition of the A matrix. This is a class method which creates and returns a ModalMvarMetrics instance containing the decomposition parameters with methods computing various connectivity metrics.

This class computes the transfer function H using modal parameters.

Currently only implemented for a single realisation (A.ndim == 3)

Parameters:

• model (sails LinearModel) – SAILS class instance containing a fitted linear model.
• sample_rate (float) – Sample rate of the dataset
• freq_vect (ndarray) – Vector of specifying which frequencies to compute connectivity at.
• sum_modes (bool) – Flag indicating whether to sum across modes (Default value = True)

Returns:

Return type:

sails.ModalMvarMetrics instance

classmethod initialise_from_modes(model, modes, sample_rate, freq_vect, mode_inds=None, sum_modes=True)

Compute some basic values from the pole-residue decomposition of the A matrix. This is a class method which creates and returns a sails.ModalMvarMetrics instance containing the decomposition parameters with methods computing various connectivity metrics.

This class computes the transfer function H using modal parameters. A sub-set of modes can be used to compute a reduced H.

Currently only implemented for a single realisation (A.ndim == 3)

Parameters:

• model (sails LinearModel) – SAILS class instance containing a fitted linear model.
• modes (sails MvarModalDecomposition) – SAILS class instance containing a fitted modal decomposition
• sample_rate (float) – The sample rate of the fitted data
• freq_vect (ndarray) – Vector of specifying which frequencies to compute connectivity at.
• mode_inds (ndarray) – Which modes to use when constructing H (Default value = None)
• sum_modes (Boolean) – Flag indicating whether to sum across modes when constructing H (Default value = True)

Returns:

Return type:

sails.ModalMvarMetrics instance

sails.mvar_metrics.modal_transfer_function(evals, evecl, evecr, nchannels, sample_rate=None, freq_vect=None)

Compute the transfer function in pole-residue form, splitting the system into modes with separate transfer functions. The full system transfer function is then a linear sum of each modal transfer function.

Parameters:

• evals (ndarray) – Complex valued eigenvalues from an eigenvalue decomposition of an MVAR parameter matrix.
• evecl (ndarray) – Complex valued left-eigenvectorsfrom an eigenvalue decomposition of an MVAR parameter matrix.
• evecr (ndarry) – Complex valued right-eigenvectorsfrom an eigenvalue decomposition of an MVAR parameter matrix.
• nchannels (int) – Number of channels in the decomposed system
• sample_rate (float) – The sampling rate of the decomposed system (Default value = None)
• freq_vect (ndarray) – Vector of frequencies at which to evaluate the transfer function (Default value = None)

Returns:

Modal transfer-function (H) of size [nchannels x nchannels x nfrequencies x nmodes]

Return type:

ndarray

Metric Mathemetical Functions

sails.mvar_metrics.sdf_spectrum(A, sigma, sample_rate, freq_vect)

Estimate of the Spectral Density as found on wikipedia and [Quirk1983].

This assumes that the spectral representation of A is invertable

Parameters:

• A (ndarray) – Matrix of autoregressive parameters, of size [nchannels x nchannels x model order]
• sigma (ndarray) – Residual covariance matrix of the modelled system
• sample_rate (float) – The samplingfrequency of the modelled system
• freq_vect (ndarray) – Vector of frequencies at which to evaluate the spectrum

Returns:

Frequenecy transform of input A, of size [nchannels x nchannels x nfrequencies]

Return type:

ndarray

sails.mvar_metrics.psd_spectrum(A, sigma, sample_rate, freq_vect)

Estimate the PSD representation of a set of MVAR coefficients as stated in [Penny2009] section 7.4.4.

This assumes that the spectral representation of A is invertable

WARNING: does not behave as expected for some data, use with caution. ar_spectrum is recommended.

Parameters:

• A (ndarray) – Matrix of autoregressive parameters, of size [nchannels x nchannels x model order]
• sigma (ndarray) – Residual covariance matrix of the modelled system
• sample_rate (float) – The samplingfrequency of the modelled system
• freq_vect (ndarray) – Vector of frequencies at which to evaluate the spectrum

Returns:

Frequenecy transform of input A, of size [nchannels x nchannels x nfrequencies]

Return type:

ndarray

sails.mvar_metrics.ar_spectrum(A, sigma, sample_rate, freq_vect)

Estimate the spectral representation of a set of MVAR coefficients as suggested by [Baccala2001], the equation without a number just below equation 13.

Parameters:

• A (ndarray) – Matrix of autoregressive parameters, of size [nchannels x nchannels x model order]
• sigma (ndarray) – Residual covariance matrix of the modelled system
• sample_rate (float) – The samplingfrequency of the modelled system
• freq_vect (ndarray) – Vector of frequencies at which to evaluate the spectrum

Returns:

Frequenecy transform of input A, of size [nchannels x nchannels x nfrequencies]

Return type:

ndarray

sails.mvar_metrics.transfer_function(Af)

Function for computing the transfer function of a system from the frequency transform of the autoregressive parameters.

Parameters:Af – Frequency domain version of parameters (can be calculated using ar_spectrum function) [nchannels x nchannels x nfreqs x nepochs] Returns:System transfer function, of size [nchannels x nchannels x nfreqs x nepochs] Return type:ndarray

Note

This function computes this the transfer function \(H\) as defined by:

\[H(f) = ( I - Af(f) )^{-1}\]

where \(Af\) is the frequency transformed autoregessive parameter matrix and \(I\) is identity.

sails.mvar_metrics.spectral_matrix(H, noise_cov)

Function for computing the spectral matrix, the matrix of spectra and cross-spectra.

Parameters:

• H (ndarray) – The transfer matrix [nchannels x nchannels x nfreqs x nepochs]
• noise_cov (ndarray) – The noise covariance matrix of the system [nchannels x nchannels x nepochs]

Returns:

System spectral matrix, of size [nchannels x nchannels x nfreqs x nepochs]

Return type:

ndarray

Note

This function computes the system spectral matrix as defined by:

\[S(f) = H(f)\Sigma H(f)^H\]

where \(H\) is the transfer function, \(\Sigma\) is the residual covariance and \(^H\) the Hermitian transpose.

sails.mvar_metrics.coherency(S)

Method for computing the Coherency. This is the complex form of coherence, from which metrics such as magnitude squared coherence can be derived

Parameters:S (ndarray) – The system spectral matrix in 3D or 4D form, of size [nchannels x nchannels x nfreqs x nepochs] Returns:Complex-valued coherency matrix, of size [nchannels x nchannels x nfreqs x nepochs] Return type:ndarray

Note

This function computes the coherency as defined by:

\[Coh_{ij}(f) = \frac{S_{ij}(f)}{ \sqrt{ | S_{ii}(f)S_{jj}(f) | } }\]

where \(S\) is the system spectral matrix.

sails.mvar_metrics.partial_coherence(inv_S)

Method for computing the Partial Coherence.

Parameters:inv_S (ndarray) – The system spectral matrix in 3D or 4D form, of size [nchannels x nchannels x nfreqs x nepochs] Returns:Complex-valued coherency matrix, of size [nchannels x nchannels x nfreqs x nepochs] Return type:ndarray

Note

This function computes the partial coherence using the inverse spectral matrix \(P = S^{-1}\) where \(S\) is the system spectral matrix. The partial coherence is then computed as:

\[PCoh_{ij}(f) = \frac{ |P_{ij}(f)|^2 }{ \sqrt{ | P_{ii}(f)P_{jj}(f) | } }\]

sails.mvar_metrics.partial_directed_coherence(Af)

Function to estimate the partial directed coherence from a set of multivariate parameters as defined in [Baccala2001].

Parameters:Af (ndarray) – Frequency domain version of parameters (can be calculated using ar_spectrum function), of size [nchannels x nchannels x nfrequencies x nrealisations] Returns:The partial directed coherence of the system. Return type:ndarray

Note

This function computes the partial coherence using the frequency transform of the autoregressive parameters \(Af\) subtracted from Identity.

\[\bar{Af}(f) = I - Af(f)\]

The partial directed coherence is then

\[PDC_{ij}(f) = \frac{ | \bar{Af}_{ij}(f)| }{ \sqrt{ \sum_i | \bar{Af}_{ij}(f) |^2 } }\]

sails.mvar_metrics.isolated_effective_coherence(Af, noise_cov)

Function for estimating the Isolated Effective Coherence as defined in [Pascual-Marqui2014].

Parameters:

• Af (ndarray) – Frequency domain version of parameters (can be calculated using ar_spectrum function) [nchannels x nchannels x nfreqs x nepochs]
• noise_cov (ndarray) – The noise covariance matrix of the system [nchannels x nchannels x nepochs]

Returns:

The isolated effective coherence.

Return type:

ndarray

sails.mvar_metrics.directed_transfer_function(H)

Method for computing the Directed Transfer Function as defined in [Kaminski1991].

Parameters:H (ndarray) – The transfer matrix [nchannels x nchannels x nfreqs x nepochs] Returns:The directed transfer function of the system. Return type:ndarray

Note

This function computes the directed transfer function from the system transfer function \(H\).

\[DTF_{ij}(f) = \frac{ | H_{ij}(f) |^2 }{ \sqrt{ \sum_j | \bar{H}_{ij}(f) |^2 } }\]

sails.mvar_metrics.geweke_granger_causality(S, H, sigma)

This function computes the Geweke-Granger causality as defined in [Barrett2010]

Parameters:

• S (ndarray) – Spectral matrix [nchannels x nchannels x nfreqs x nepochs]
• H (ndarray) – The transfer matrix [nchannels x nchannels x nfreqs x nepochs]
• sigma (ndarray) – Residual noise covariance matrix

Returns:

The Geweke-Granger causality

Return type:

ndarray

Plotting

sails.plotting.root_plot(rts, ax=None, figargs={}, plotargs={})

Plot a set of roots (complex numbers).

Parameters:

• rts (ndarray_like) – Roots to plot
• ax (matplotlib axes handle) – Optional Axes on which to place plot. (Default value = None)
• figargs (dict) – extra arguments to pass to plt.figure (Default value = dict())
• plotargs (dict) – extra arguments to pass to plt.plot (Default value = dict())

Returns:

Axes object on which plot was drawn

Return type:

Axes

sails.plotting.plot_vector(metric, x_vect, y_vect=None, x_label=None, y_label=None, title=None, labels=None, line_labels=None, F=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, triangle=None, diag=False, thresh=None, two_tail_thresh=False, font_size=10, use_latex=False)

Function for plotting frequency domain connectivity at a single time-point.

Parameters:

• metric (ndarray) – matrix containing connectivity values [nsignals x signals x frequencies x participants] in which the first dimension refers to source nodes and the second dimension refers to target nodes
• x_vect (ndarray) – vector of frequencies to label the x axis
• y_vect (ndarrat) – vector containing the values for the y-axis
• x_label (string [optional]) – label for the x axis (Default value = None)
• y_label (string [optional]) – label for the y axis (Default value = None)
• title (string [optional]) – title for the figure (Default value = None)
• labels (list) – list of node labels for columns and vectors (Default value = None)
• line_labels (list) – list of labels for each separate line (participant dimension in metric) (Default value = None)
• F (figurehandle [optional]) – handle of existing figure to plot within (Default value = None)
• triangle (string [optional]) – string to indicate whether only the ‘upper’ or ‘lower’ triangle of the matrix should be plotted (Default value = None)
• diag (bool [optional]) – flag to indicate whether the diagonal elements should be plotted (Default value = False)
• thresh (ndarray [optional]) – matrix containing thresholds to be plotted alongside connectivity values [nsignals x nsignals x frequencies] (Default value = None)
• two_tailed_thresh (bool [optional]) – flag to indicate whether both signs (+/-) of the threshold should be plotted
• font_size (int [optional]) – override the default font size
• use_latex (bool) – Flag to indicate whether to render text in latex (Default value = False)

Returns:

Figure handle containing the plot

Return type:

matplotlib figure handle

sails.plotting.plot_matrix(metric, x_vect, y_vect, x_label=None, y_label=None, z_vect=None, title=None, labels=None, F=None, vlines=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, font_size=8, use_latex=False, diag=True)

Function for plotting frequency domain connectivity over many time points

Parameters:

• metric (ndarray) – matrix containing connectivity values [nsignals x signals x frequencies x participants] in which the first dimension refers to source nodes and the second dimension refers to target nodes
• x_vect (1d array) – vector of frequencies to label the x axis
• y_vect (1d array [optional]) – vector containing the values for the y-axis
• z_vect (1d array [optional]) – vector containing values for the colour scale (Default value = None)
• x_label (string [optional]) – label for the x axis (Default value = None)
• y_label (string [optional]) – label for the y axis (Default value = None)
• title (string [optional]) – title for the figure (Default value = None)
• labels (list) – list of node labels for columns and vectors (Default value = None)
• F (figurehandle [optional]) – handle of existing figure to plot within (Default value = None)
• vlines (list) – List of x-axis values to plot a dashed vertical line (Default value = None)
• cmap (matplotlib colormap [optional]) – matplotlib.cm.<colormapname> to use for colourscale (redundant for plot vector??) (Default value = plt.cm.jet)
• font_size (int [optional]) – override the default font size
• use_latex (bool) – Flag indicating whether to render text in latex (Default value = False)
• diag (bool) – Flag indicating whether to plot the diagonal subplots (Default value = True)

Returns:

Figure handle containing the plot

Return type:

matplotlib figure handle

class sails.circular_plots.CircosHandler

Handler for producing plots using the Circos tool.

Note that if you are creating plots for publication using Circos, you should cite the relevant publication: [Krzywinski2009].

class sails.circular_plots.CircosGroup(name, c_idx=0, n_idx=0, circoscolour='unknown', rgbcolour=None, netmatreverse=False)

Support routines

Abstract Classes

class sails.modelfit.AbstractLinearModel

This is a base class defining data storage and generic methods for LinearModel classes. New LinearModels should inherit from AbstractLinearModel and overload the fit_model method

companion

Parameter matrix in companion form

compute_diagnostics(data, compute_pc=False)

Compute several diagnostic metrics from a fitted MVAR model and a dataset.

Parameters:

• data (ndarrat) – The data to compute the model diagnostics with, typically the same data used during model fit.
• compute_pc (bool) – Flag indicating whether the compute the percent consistency, this is typically the the longest computation of all the metrics here (Default value = False)

Returns:

object containing LL, AIC, BIC, stability ratio, durbin-watson, R squared and percent consistency measures

Return type:

sails.ModelDiagnostics

get_residuals(data, mode='valid')

Returns the prediction error from a fitted model. This is a wrapper function for get_residuals()

Parameters:data (ndarray) – Data to compute the residuals of, of size [nchannels x nsamples x nrealisations] Returns:Residual data Return type:ndarray

nsignals

Number of signals in fitted model

order

Order of fitted model

class sails.mvar_metrics.AbstractMVARMetrics

This is a base class for computing various connectivity metrics from fitted autoregressive models

PSD

Power Spectral Density matrix

S

Spectral Matrix of shape [nsignals, nsignals, order, epochs]

coherency

Complex-valued coherency

d_directed_transfer_function

Direct directed transfer function

directed_transfer_function

Directed transfer function

ff_directed_transfer_function

Full-frequency directed transfer function

get_node_weight(metric='S')

Compute the node weight for a given connectivity metric.

Parameters:metric (str) – Name of the connectivity metric to plot (Default value = ‘S’) Returns:Node-weight matrix of size [nchannels x nchannels] Return type:ndarray

get_spectral_generalisation(metric='S')

Compute the spectral generalisation matrix for a given metric. Returns the spatial similarity of different frequencies.

Parameters:metric (str) – Name of the connectivity metric to plot (Default value = ‘S’) Returns:Matrix of spectral generalisation [nfrequencies x nfrequencies] Return type:ndarray

geweke_granger_causality

Geweke granger causality

imaginary_coherence

Imaginary part of complex coherency

inv_S

Inverse Power Spectrl Density matrix

isolated_effective_coherence

Isolated effective coherence

magnitude_squared_coherence

Magnitude squared coherence

partial_coherence

Partialled magnitude-squared coherence

partial_directed_coherence

Partial directed coherence

phase_coherence

Phase angle of complex coherency

plot_diags(metric='S', inds=None, F=None, ax=None)

Helper function for plotting the diagonal part of a connectivty metric. Calls sails.plotting.plot_diagonal.

Parameters:

• metric (str) – Name of the connectivity metric to plot (Default value = ‘S’)
• F (matplotlib figure handle) – Handle for matplotlib figure to plot in (Default value = None)
• ax (matplotlib axes handle) – Handle for matplotlib axes to plot in (Default value = None)

plot_summary(metric='S', ind=1)

Helper function for plotting a summary of a connectivty metric. Calls sails.plotting.plot_summary.

Parameters:

• metric (str) – Name of the connectivity metric to plot (Default value = ‘S’)
• ind (int) – Index into frequency dimension to plot connectivity matrix (Default value = 1)

Statistical Functions

sails.stats.find_cov_multitrial(a, b, ddof=1)

Estimate the average covariance between two datasets across many trials

Parameters:

• a (ndarray) – a multivariate time series of shape nsignals x nsamples x ntrials
• b (ndarray) – a multivariate time series of shape nsignals x nsamples x ntrials
• ddof (int) – if None, no normalisation will be performed. If an integer, covariance will be normalised by (nsamples - ddof) (Default value = 1)

Returns:

A covariance matrix between a and b of shape nsignals x nsignals

Return type:

type

sails.stats.find_cov(a, b, ddof=1)

Estimate the covariance between two datasets.

Parameters:

• a (ndarray) – a multivariate time series of shape nsignals x nsamples
• b (ndarray) – a multivariate time series of shape nsignals x nsamples
• ddof (int) – if None, no normalisation will be performed. If an integer, covariance will be normalised by (nsamples - ddof) (Default value = 1)

Returns:

A covariance matrix between a and b of shape nsignals x nsignals

Return type:

type

sails.stats.mutual_information(A, B, r=None, nbins=21, base=<ufunc 'log2'>)

Calculate the mutual information between two signals.

Parameters:

• A (ndarray) – First signal, of shape [channels x samples]
• B (ndarray) – Second signal, of shape [channels x samples]
• r (tuple (min,max)) – Range of the random variables used for bin calculations. If not set, takes the minimum and maximum values from the combination of A and B (Default value = None)
• nbins (int) – Number of bins to compute the frequencies within (Default value = 21)
• base (function) – Base to use for MI calculation - defaults to np.log2

Returns:

the mutual information within each signal passed

Return type:

type

sails.stats.durbin_watson(residuals, step=1)

Test for autocorrelation in the residuals of a model (see [Durbin1950] and [Durbin1951])

The result is a value between 0 and 4. A value of 0 indicates a positive autocorrelaiton, 2 indicates no correlation and a value of 4 indicates a negative autocorrelation.

If many trials are found the statistic is estimated for each trial separately.

Parameters:

• residuals (ndarray) – Residuals of shape [channels, samples, ntrials] If only two dimensions, will be treated as [channels, samples] If only one dimension, will be treated as [samples]
• step (int) – step-size to use through residuals when calculating DW. Defaults to 1.

Returns:

durbin-watson index

Return type:

type

sails.stats.percent_consistency(data, residuals)

Estimate how well a model retains the correlational structure of a dataset.

This uses the percent consistency measure outlined in [Ding2000]. All auto and cross correlations within a multivariate data set are computed and compared to the auto and cross correlations of the data as predicted by a model (in this case the data - model residuals, as these are easily obtained).

Parameters:

• data (ndarray) – An array containing the observed data, of shape [nsignals x nsamples x ntrials]
• residuals (ndarray) – An array containing the remaining variance of data after model fitting, of shape [nsignals x nsamples x ntrials]

Returns:

A value between 0 and 100 indicating the percentage of the autocorrelation in the data that is retained in the model.

Return type:

type

sails.stats.rsquared_from_residuals(data, residuals, per_signal=False)

Estimate the variance explained by a model from the original data and the model residuals

If data and residuals are 2D, they will be treated as [channels, nsamples]

If data and residuals are 1D, they will be treated as [nsamples]

Parameters:

• data (ndarray) – Array containing the original data, of shape [channels, nsamples, ntrials]
• residuals (ndarray) – Array containing the model residuals, of shape [channels, nsamples, ntrials]
• per_signal (bool) – Boolean indicating whether the variance explained should be computed separately per channel (Default value = False)

Returns:

Value between 0 and 1 indicating the proportion of variance explained by the model.

Return type:

type

sails.stats.mahalanobis(X, Y=None, sigma=None)

Estimate the Mahalanobis distance of each point in an array across the samples dimension.

Parameters:

• X (ndarray) – Array containing data of points: shape [channels x samples]
• Y (ndarray (Optional)) – Optional Array containing data of second points: shape [channels x samples]. If not provided, the origin will be used as the point for comparison. (Default value = None)
• sigma (ndarray (Optional)) – A precomputed covariance matrix to use in the calculation, of shape [channels x channels]. If not provided, sigma will be calculated from X (Default value = None)

Returns:

The Mahalanobis distance for each sample, of shape [samples]

Return type:

type

sails.stats.profile_likelihood(X)

Find the ‘elbow’ point in a curve.

This function uses the automatic dimensionality selection method from [Zhu2006] to estimate the elbow or inflection point in a given curve.

Parameters:X (ndarray) – 1d array containing the curve to be analysed, of shape [samples] Returns:The log likelihood that each point on the curve separates the distributions to the left and right, of shape [samples] Return type:type

sails.stats.stability_index(A)

Compute the stability index as defined in [Lutkephol2006] pages 15-16.

This is a proxy for the stationarity assumption of MVAR models, if the magnitude of the largest principal component of the model parameters is less than 1, the model is stable and thus, stationary.

Parameters:A (ndarray) – The MVAR parameter matrix of size [nchannels x nchannels x model_order] Returns:The stability index Return type:type

sails.stats.stability_ratio(A)

Compute the stability ratio as a measure of stationarity

A stronger test for stationarity than the stability index. Computes the ratio between the largest eigenvalue of the parameters and all the others. If the largest parameter is larger than all other parameters (indicated by a value < 1) we have a super-stable system.

Parameters:A (ndarray) – The MVAR parameter matrix of size [nchannels x nchannels x model_order] Returns:The stability ratio Return type:type

sails.modelvalidation.approx_model_selection_criteria(A, EF, method='ss')

Estimate Akaike’s Information Criterion and Bayes Information Criterion from a model fit using an approximation of the loglikelihood.

Two approximations are available:

‘ss’: The sum square of the residuals. This approximation is outlined in [Burnham2002] and [Chatfield2003]. ‘det’: The determinant of the residual covariance matrix.

The sum squares method is used as default.

Parameters:

• A (ndarray) – Array containing the parameters of a fitted model, of shape [channels x channels x maxorder]
• EF (ndarray) – Array containing the residuals of the fitted model, of shape [channels x samples x trials]
• method ({'ss','det'}) – string indicating which method should be used to approximate the loglikelihood (‘ss’ or ‘det’) (Default value = ‘ss’)

Returns:

tuple of (aLL, AIC, BIC) aLL: Approximation of the loglikelihood (determined by the method defined above) AIC: Estimate of Akaike’s Information Criterion BIC: Estimate of the Bayesian Information Criterion

Return type:

tuple

Data Processing Functions

sails.modal.adjust_phase(x)

Method for normalising eigenvector matrices to meet the conditions in eqn 21 of [Neumaier2001].

Parameters:x (ndarray) – Vector of complex numbers to adjust the phase of Returns:Vector of numbers with adjusted phase Return type:type

sails.orthogonalise.symmetric_orthonormal(A, maintain_mag=False)

Implement the [Colclough2015] algorithm for multivariate leakage correction in MEG. Also see the OHBA toolbox (https://github.com/OHBA-analysis/MEG-ROI-nets) and http://empslocal.ex.ac.uk/people/staff/reverson/uploads/Site/procrustes.pdf

Parameters:

• A (ndarray) – Data matrix to orthogonalise, of size [channels, samples]
• verbose (bool) – Flag indicating whether to show detailed output (Default value = False)
• maintain_mag (bool) – Flag indicating whether to maintain magnitudes (Default value = False)

Returns:

• data (numpy.ndarray) – [channels, samples] of orthogonalised data
• W (numpy.ndarray)

sails.orthogonalise.innovations(A, order, mvar_class=None, verbose=False)

Implementation of the innovations-orthogonalisation defined in [Pascual-Marqui2017] . This creates the mixing matrix on the input data after variance explained by an MVAR model is removed. The orthogonalisation is defined on the residuals but applied to the raw data.

Parameters:

• A (ndarray) – Data matrix to orthogonalise, of size [channels, samples]
• order (int) – Model order used during MVAR model fit
• mvar_class (sails LinearModel class) – Model class to do model fitting (Default value = None)
• verbose (bool) – Flag indicating whether to show detailed output (Default value = False)

Returns:

Orthogonalised data matrix, of size [channels, samples]

Return type:

ndarray

Tutorial helper functions

sails.tutorial_utils.generate_pink_roots(power, order=24)

Generates the poles of a system with a pink spectrum according to eqn 116 in [Kasdin1995]

Parameters:

• power (float) – Power of pole must be between -2 and 2.
• order – Order of system to generate (Default value = 24)

Returns:

Roots of the system (length of order)

Return type:

ndarray

sails.tutorial_utils.model_from_roots(roots)

Sets up a univariate model based on the roots of a polynomial

Parameters:roots (ndarray) – Vector of polynomial roots. Returns:Model fit class containing the parameters defined by the input roots Return type:AbstractLinearModel

sails.tutorial_utils.find_support_path()

Returns the path to the support files directory

sails.tutorial_utils.find_example_path()

Returns the path to the tutorial example data directory (if available)

sails.tutorial_utils.set_sail()

References

[Baccala2001]

Baccala L.A. and Sameshima K. Partial directed coherence: a new concept in neural structure determination. Biological Cybernetics. 84. 463-474. https://doi.org/10.1007/PL00007990

[Barrett2010]

Barrett A.B., Barnett L. and Seth A.K. Multivariate Granger causality and generalized variance. Physical Review E. E81, 041907. http://dx.doi.org/10.1103/PhysRevE.81.041907

[Burnham2002]

Burnham K.P. and Anderson D.R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer. ISBN: 9780387224565

[Chatfield2003]

Chatfield C. The Analysis of Time Series: An Introduction. Chapman and Hall. ISBN: 9780429208706

[Colclough2015]

Colclough G.L., Brookes M.J., Smith S.M., Woolrich M.W. A symmetric multivariate leakage correction for MEG connectomes. NeuroImage. 117. 439-448. http://dx.doi.org/10.1016/j.neuroimage.2015.03.071

[Ding2000]

Ding M., Bressler S.L., Weiming Y. and Liang H. Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variability assessment. Biological Cybernetics. 83(1) 35-45. https://doi.org/10.1007/s004229900137

[Durbin1950]

Durbin, J., & Watson, G. S. (1950). Testing For Serial Correlation In Least Squares Regression. I. Biometrika, 37(3–4), 409–428. https://doi.org/10.1093/biomet/37.3-4.409

[Durbin1951]

Durbin, J., & Watson, G. S. (1951). Testing For Serial Correlation In Least Squares Regression. II. Biometrika, 38(1–2), 159–178. https://doi.org/10.1093/biomet/38.1-2.159

[Fasoula2013]

Fasoula A, Attal Y and Schwartz D. Comparative performance evaluation of data-driven causality measures applied to brain networks. Journal of Neuroscience Methods. 215(2). 170-189. https://doi.org/10.1016/j.jneumeth.2013.02.021

[Kaminski1991]

Kaminski M.J. and Blinowska K.J. A new method of the description of the information flow in the brain structures. Biological Cybernetics. 65(3) 203-210. http://dx.doi.org/10.1007/BF00198091

[Kasdin1995]

Kasdin, N.J. Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^alpha Power Law Noise Generation, Proceedings of the IEEE, Vol. 83, No. 5, 1995, pp. 802-827.

[Korzeniewska2003]

Korzeniewska A, Mańczakb M, Kamiński M, Blinowska K.J. and Kasicki S. Determination of information flow direction among brain structures by a modified directed transfer function (dDTF) method. Journal of Neuroscience Methods. 125(1-2). 195-207. https://doi.org/10.1016/S0165-0270(03)00052-9

[Krzywinski2009]

Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, Jones S.J. and Marra M.A. Circos: an Information Aesthetic for Comparative Genomics. Genome Res. 19(9). 1639-1645. http://dx.doi.org/10.1101/gr.092759.109

[Lutkephol2006]

Lutkepohl H. New Introduction to Multiple Time Series Analysis. Springer. ISBN: 9783540262398

[Marple1987]

Marple Jr, L.S. Digital Spectral Analysis: With Applications. Prentice-Hall. ISBN: 9780132141499

[Neumaier2001]

Neumaier A. and Schneider T. Estimation of parameters and eigenmodes of multivariate autoregressive models. ACM Transactions on Mathematical Software. 27(1). 27-57. http://dx.doi.org/10.1145/382043.382304

[Pascual-Marqui2014]

Pascual-Marqui R.D., Biscay R.J., Bosch-Bayard J., Lehmann D., Kochi K., Kinoshita T., Yamada N. and Sadato N. Assessing direct paths of intracortical causal information flow of oscillatory activity with the isolated effective coherence (iCOH). Frontiers in Human Neuroscience. 8. 448. http://dx.doi.org/10.3389/fnhum.2014.00448

[Pascual-Marqui2017]

Pascual-Marqui R.D., Biscay R.J., Bosch-Bayard J., Faber P., Kinoshita T., Kochi K., Milz P., Nishida K. and Yoshimura M. Innovations orthogonalization: a solution to the major pitfalls of EEG/MEG “leakage correction”. Biorxiv. https://doi.org/10.1101/178657

[Penny2009]

Penny W.D. Signal Processing Course. https://www.fil.ion.ucl.ac.uk/~wpenny/course/course.html

[Quinn2020]

Quinn A.J. and Hymers M. SAILS: Spectral Analysis In Linear Systems. Journal of Open Source Software, 5(47), 1982. https://doi.org/10.21105/joss.01982

[Quinn2021]

Quinn A.J., Green G.G.R. and Hymers, M. Delineating between-subject heterogeneity in alpha networks with Spatio-Spectral Eigenmodes. NeuroImage, 118330. https://doi.org/10.1016/j.neuroimage.2021.118330

[Schelter2006]

Schelter B, Winterhalder M, Eichler M, Peifer M, Hellwig B, Guschlbauer B, Hermann-Lucking C, Dahlhaus R and Timmer J. Testing for directed influences among neural signals using partial directed coherence. Journal of Neuroscience Methods. 152(1-2). 210-219. https://doi.org/10.1016/j.jneumeth.2005.09.001

[Quirk1983]

Quirk M. and Liu B. Improving resolution for autoregressive spectral estimation by decimation. IEEE Transactions on Acoustics, Speech and Signal Processing. 31(3) 630-637. http://dx.doi.org/10.1109/TASSP.1983.1164124

[Zhu2006]

Zhu M and Ghodsi A. Automatic dimensionality selection from the scree plot via the use of profile likelihood. Computational Statistics & Data Analysis. 51(2) 918-930. https://doi.org/10.1016/j.csda.2005.09.010

Citing SAILS

If this package has contributed to your work, please include the following citations:

Quinn, A. J., & Hymers, M (2020). SAILS: Spectral Analysis In Linear Systems. Journal of Open Source Software, 5(47), 1982,
https://doi.org/10.21105/joss.01982

Quinn, A. J., Green, G. G. R., & Hymers, M. (2021). Delineating between-subject heterogeneity in alpha networks with Spatio-Spectral Eigenmodes. NeuroImage, 118330.
https://doi.org/10.1016/j.neuroimage.2021.118330

@article{Quinn2020,
  doi = {10.21105/joss.01982},
  url = {https://doi.org/10.21105/joss.01982},
  year = {2020},
  publisher = {The Open Journal},
  volume = {5},
  number = {47},
  pages = {1982},
  author = {Andrew J. Quinn and Mark Hymers},
  title = {SAILS: Spectral Analysis In Linear Systems},
  journal = {Journal of Open Source Software}
}

@article{Quinn2021,
  doi = {10.1016/j.neuroimage.2021.118330},
  url = {https://doi.org/10.1016/j.neuroimage.2021.118330},
  year = {2021},
  month = jul,
  publisher = {Elsevier {BV}},
  pages = {118330},
  author = {Andrew J. Quinn and Gary G.R. Green and Mark Hymers},
  title = {Delineating between-subject heterogeneity in alpha networks with Spatio-Spectral Eigenmodes},
  journal = {{NeuroImage}}
}