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
    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.