Chapter 5.4. Model Fitting and Validation

Once the data has been pre-processed, we can proceed to Model Fitting and Validation. SIFT currently supports parametric VAR modeling using the Vieira-Morf (Marple, 1987), or ARFIT algorithm (Schneider and Neumaier, 2001). ARFIT must be downloaded separately and installed in /external/arfit/). SIFT currently supports time-varying parametric VAR modeling either through the use of sliding-window adaptive VAR modeling (est_fitMVAR()) or recursive-least-squares Kalman filtering (est_fitMVARKalman()).

Model Fitting and Validation can also be started from the command line:

EEG = pop_est_fitMVAR(EEG);

You should now see a GUI similar to that displayed in the figure below. Here we can select the MVAR algorithm, choose the sliding window length and step size, and specify the model order. ARFIT uses a modified least-squares approach, while Vieira-Morf uses a multichannel geometric-mean non-least-squares lattice approach to solve for the model coefficients. ARFIT includes an additional term for the process mean, whereas Vieira-Morf assumes a zero-mean process. The current implementation of ARFIT is faster than that of the Vieira-Morf algorithm, although Vieira-Morf returns slightly better coefficient estimates. For a detailed performance comparison between these and other estimators, see Schlögl (2006). In this example, we will use the Vieira-Morf algorithm. We will also use a window length of 0.4 sec (400 ms) with a step size of 0.01 sec (10 ms). This is approximately 1 cycle of a 2.5 Hz oscillation and should allow us to adequately model information flow from the lowest frequency band of interest (delta) up to our Nyquist frequency of 128 Hz. Consult the theory section below for more details on selecting an appropriate window length.

Here we can choose to calculate one or more of the information criteria listed in the section 3.5. over a range of model orders (p_min to p_max). If more than one window is available, the information criteria are calculated for each window, and histograms will be generated showing the distribution of optimal model orders for all windows. If we have a large number of windows and are pressed for time, we can specify a random percentage of windows for which to calculate the criteria (% windows to sample). Bear in mind, however, that increasing the number of windows used will result in a better estimate of the empirical distribution of the information criteria across windows. Additionally, for a fast, approximate estimate of the information criteria, we can choose to downdate the model (Neumaier and Schneider, 2001). Rather than fitting p_max – p_min VAR models of increasing order, this fits a single VAR[p_max] model, and downdates the noise covariance matrix to obtain approximate estimates of the information criteria.

Your GUI should appear as shown in the figure below with options set as in the table below:

Option	Value
Select MVAR Algorithm	Vieira-Morf
Window length (sec)	0.4
Step Size (sec)	0.01
Model Order	[1 30]
Downdate	True
Model Order	[1 30]
Windows to sample	100

Figure caption. AMVAR model fitting GUI generated by pop_est_fitMVAR(). We have selected a window length of 0.4 sec, a step size of 0.01 second, and specified a model order range (for order selection) of 1 to 30.

Before we move forward, let's see how to select the optional window length.

5.4.1. Theory: selecting a window length

There are several important considerations that can aid us in choosing an appropriate window length. An appropriate window length often requires a compromise between one or more of the following considerations:

1. Local Stationarity
2. Temporal Smoothing
3. Sufficient amount of data
4. Process Dynamics and Neurophysiology

5.4.1.1. Local Stationarity

In section 3.4., we discussed issues surrounding non-stationary EEG data and introduced the concept of using a sliding window to fit VAR models to locally stationary data. It is thus important that our window be small enough to ensure local stationarity. As we shall see in later in this page, validating our model and plotting the stability index for a chosen window size can give us a sense as to whether our window is sufficiently small to ensure local stationarity.

5.4.1.2. Temporal Smoothing

A sliding-window analysis is inherently a temporal smoothing operation. Larger windows result in model coefficients being aggregated over increasingly longer segments of data and, therefore, result in increased temporal smoothing. If there are rapid transient changes in process dynamics, choosing a too-large window may obscure the fine temporal structure of information flow. When multiple trials are available, a useful heuristic proposed by Ding et al. (2000b) for obtaining a rough upper limit on the window length is to plot the bivariate correlation for a few pairs of variables across all trials, beginning with a 1-point window and successively averaging correlation within trials and across the window for increasingly larger windows. An illustration from Ding et al. (2000b) is reproduced in the figure below. Note that with the 1-point window, there are large fluctuations in correlation structure (covariance non-stationarity). As we increase the window length, we get increased temporal smoothing. In this case, a reasonable window length might be 10 points, since it reduces local covariance non-stationarity (local fluctuations in cross-correlation) while still preserving some of the temporal dynamics of interest (namely the changes in correlation structure). Of course, this suggested window length is completely application-specific; one should select a window tailored to their specific data.

Figure caption. Cross-correlation between two intracranial EEG time series averaged over increasing window lengths (1 point, 10 points, 20 points) and plotted as a function of time. Figure reproduced from Ding et al. (2000b).

5.4.1.3. Sufficient amount of data

In section 3.4.1. we noted that a minimum of M²p data points are required to fit an M-dimensional VAR[p] model. We also stated that, in practice, we would like to have ten times that many data points to ensure an unbiased fit (if the model order is 15, we want at least 150 data points to fit the model). This leads us to the rule-of-thumb equation:

or, equivalently,

where W is the window length in points and N is the total number of trials available. SIFT performs checks on parameters (est_checkMVARParams()) and will let you know if the selected window length is sub-optimal (as well as recommend a better window length).

5.4.1.4. Process dynamics and neurophysiology

In section 3.5., we discussed how, when selecting an appropriate model order, one should take into account the temporal dynamics and neurophysiology of the data. The same principles hold for window length selection. Although, with a large number of trials, we could theoretically fit our model using a window length as short at p+1 sample points long, we must consider that all interactions being modeled must occur within the selected window. In general, if we expect a maximum interaction lag of seconds between any two brain processes, we should make sure to select a window length of W .

Furthermore, if we are interested in frequency-domain quantities, we should consider the time-frequency uncertainty principle, which states that every increase in temporal resolution leads to a concomitant decrease in frequency resolution. In general, a useful rule-of-thumb is to ensure that the window is long enough to span approximately one cycle of the lowest frequency of interest.

5.4.2. Selecting the model order

Now that we have chosen our VAR algorithm, window length, and step size, we can proceed to model order analysis. Upon pressing OK in the previous GUI. You should see a warning window pop up below. It shows the results of a sanity check that evaluates the ratio of parameters to data points, calculates the number of estimable frequencies, checks the time-frequency product, and performs other relevant checks. Information is displayed for each condition, along with suggestions on optimal parameters to use if your parameter selections are sub-optimal. Here, we are being warned that the ratio of free parameters to data points is greater than 0.1, which may cause concern. This is because our upper model order of 30 is quite large. Let’s go ahead and ignore this error (we are likely to use a much lower model order when we fit our final model). The same reasoning goes for the warning about a short window size (0.4 seconds). Simply press OK.

Figure caption. This GUI shows the results of a sanity check performed on the specified model parameters (this is always performed before model fitting).

A progress bar should now indicate our progress for each condition (RespCorr, RespWrong). When this is complete, you should see the resulting figures shown below. On the left is the result of the model order selection for RespWrong, and on the right is the result for RespCorr. The top panel shows the information criteria (averaged across windows) plotted versus the model. The vertical lines indicate each criterion's average optimal model order (model order that minimizes the information criterion). The lower array of histograms shows the distribution of optimal model orders for all windows for each criterion. Note that, as mentioned in section 3.5., for many windows, the AIC and FPE criteria do not appear to exhibit a clear minimum across the specified model range. In contrast, SBC shows a clear minimum peaking around p=5 (which is likely too low given that we will only be able to estimate 2.5 spectral peaks for each pair of variables), while HQ shows a clear minimum around p=10 for RespWrong and RespCorr. Note, however, that in both RespWrong and RespCorr, the upper limit on the model order selection criteria is approximately p=15. If we click on the top panel of RespCorr, we get an expanded view of the panel. Likewise, clicking on the histogram for HQ pops up an expanded view of the histogram. The shaded region indicates the range where 90% of the optional model order lies, and the text "20+-6" indicates the mean and standard deviation for each model order selection criteria. Note that although the minimum for the HQ criterion (purple) is at p=13, the upper limit of the “elbow” for the HQ criterion is around p=15 or p=16. It also appears that AIC/FPE indicates larger values of about 20 compared to SBC and HQ. From this, we might conclude that a suitable and safe model order for all windows and conditions is p=15.

Figure caption. Results of model order selection for RespWrong (left) and RespCorr (right). The top panel plots the information criteria versus model order and marks the average optimal model order for the selected criteria.

Figure caption. Left: Close-up view of SBC (blue), HQ (purple), FPE (yellow), AIC (red; overlapping with FPE) information criteria plotted versus model order. Note that FPE and AIC plots are almost identical. Right: Distribution over all windows of optimal model order using HQ information criterion. Vertical markers denote the distribution mean. Note the distribution is somewhat bimodal, with one peak around 9 and another around 14.

5.4.3. Fitting the final model

Then, the following GUI pops up.

Figure caption. GUI that pops up at the end of model order estimation.

We are now ready to fit the model. You may press OK on the previous query window, call menu item Tools > SIFT > Model fitting and validation > Fit AMVAR model or use the command line call below:

EEG = pop_est_fitMVAR(EEG);

We now set the model order option to 15 as per the previous discussion. The final set of parameters are shown in the table below:

Option	Value
Select MVAR algorithm	Vieira-Morf
Window length (sec)	0.4 (400 ms)
Step size (sec)	0.01 (10 ms)
Model order	15

Figure caption. This GUI allows entering parameters for model order estimation. Our final set of selected parameters for model fitting. Note that we have selected the Vieira-Morf algorithm, a window length of 0.4 seconds with a step size of 0.01 sec, and a model order of 15. Upon clicking OK, a progress bar will show us the status of the model-fitting algorithm.

Click OK to continue. SIFT sanity check should proceed and generate no warnings or errors indicating we chose a valid set of model parameters, as shown below.

Figure caption. SIFT sanity check.

If you do not wish to use the GUI, the same may be achieved from the command line using:

EEG = pop_est_fitMVAR( EEG, 'nogui', 'Algorithm', 'Vieira-Morf', 'ModelOrder', 15, 'WindowLength', 0.4, 'WindowStepSize', 0.01, 'verb', 1); 

For each condition, the VAR[15] model will now be fit for each of the 238 windows. We should now see a progress bar indicating the model fitting progress for each condition. Depending on the speed of your computer, this may take a while, so you might want to take a break or do a little yoga. If you have little computer memory or processor speed issues, you can increase the step size to 0.03 sec. This will reduce the computation time demands while still producing similar results as in the remainder of this tutorial.

5.4.4. Validating the fitted model

After you are refreshed from that Yoga session and the model has been fit for each condition, we will need to validate our fitted model. Select menu item Tools > SIFT > Model Fitting and Validation > Validate Model or type on the command line:

 pop_est_validateMVAR(EEG);

You should now be presented with the GUI shown in the figure below. Here, we have the option to check the whiteness of the residuals, percent consistency, and model stability for each (or a random subset) of our windows. As we discussed in section 3.6., residual whiteness tests include portmanteau and autocorrelation tests for correlation in the residuals of the model (which could indicate a poor model fit). Here, we have the option to choose from the Ljung-Box, Box-Pierce, and Li-McLeod multivariate portmanteau tests and a simple autocorrelation function test. Percent consistency denotes the fraction of the correlation structure of the original data that is captured by the model, while model stability performs an eigenvalue test for the stability/stationarity of the process. The options for this GUI should be set as shown in the table below:

Option	Value
Check Whiteness of Residuals	checked
significance level	0.05
Check percent consistency	checked
Check model stability	checked
% windows to sample	100

Figure caption. Model Validation GUI generated by pop_est_validateMVAR(). Here, we can choose to check the whiteness of residuals, percent consistency, and model stability for all (or some random subset) of windows. In this example, we have chosen a significance threshold of p<0.05 for our whiteness tests.

Click OK to continue. You should now see a sequence of progress bars for each condition. This may take a while.

Once the model validation routines have been completed, you should see the results shown for each of the two conditions, as in the figure below. The top panel of each figure shows the results of the whiteness tests as a function of the window index (sorted in order of temporal precedence). For the portmanteau tests, we have plotted the p-value for acceptance of the null hypothesis of correlated residuals (namely 1-p where p is the p-value for rejection of the null hypothesis). Values greater than 0.05 (blue dashed line) indicates the residuals are white at the p<0.05 level. For the ACF test (red) we have plotted the probability of an observed ACF coefficient to be within the expected interval for white residuals. Values greater than 0.95 indicate the residuals are white at the p<0.05 level. The fraction of windows that pass the whiteness test is noted in the legend. Note that the ACF tests classify all windows as having white residuals, while the portmanteau tests (which are much more conservative) indicate that most windows are white. The fact that a range of windows near the end of the epoch marginally fail the portmanteau tests may indicate that we may want to use a slightly larger model order (e.g., p=16) or perhaps a smaller window size (to improve local stationarity).

The middle panel shows the percent consistency plotted versus the increasing window index. Note that the PC is reliably high (µ≈87%), suggesting a reasonable model fit.

The lower panel shows the stability index for each window. Values above or near 0 indicate an unstable (and possibly non-stationary) model. In this case, we might try some additional preprocessing or shorten the window length to improve the local stationarity of the data. In our example, the stability index is reliably low, indicating a stable/stationary model.

The validation checks all indicate a reasonably fit model (although there may be room for improvement of the fit). Assuming we are comfortable with this we can now proceed to spectral/connectivity estimation and visualization.

Figure caption. Results of model validation for RespWrong (top) and RespCorr (bottom) conditions. For each condition, a validation statistic is plotted versus the window index (sorted in order of temporal precedence). If only one window is available, bar plots are generated instead. The top panel shows the significance level for rejecting the hypothesis of correlated residuals. For portmanteau tests (LMP, Box-Pierce, Ljung-Box), a value greater than 0.05 (dashed line) and a value greater than 0.95 indicate white residuals at the p > 0.05 level. The middle panel shows the percent consistency. The bottom panel shows the stability index.