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