Optimizing event-related potential based brain–computer interfaces: a systematic evaluation of dynamic stopping methods

All methods are conceptually described here, within the above framework. A more detailed description is given if the implementation deviates from the original implementation. For a full account, the reader is referred to the original publications. All evaluations are done in the context of a text spelling application and thus by optimizing the hyperparameter S of a stopping method with respect to the symbols per minute (SPM) metric (see section 2.4 and 2.5). The following methods are investigated:

No stopping. In this condition, no early stopping is performed. It uses all available iterations and is therefore the performance baseline. No parameters need to be estimated. Formally, this means that Θ = {J}, the maximum number of iterations possible. F(D^test, Θ, t, j) returns $c^w_{j,t}$ (according to equation (3)) when j = J.

Fixed (optimal). The simplest way of optimizing the BCI is to estimate an optimal fixed number of iterations on the calibration data and to apply this for the online data. Because the number of iterations can be set in virtually any BCI system, this method requires no implementation effort. The number of iterations j† which results in the highest SPM on the offline/training data is selected for online use. Θ is thus defined as Θ = {j†}. As j† is directly optimized it equals S. During online use, F(D^test, Θ, t, j) returns $c^w_{j,t}$ (equation (3)) when j = j†. For the fixed condition, the optimal number of iterations is estimated once from all of the training data; it serves as a evaluation of the post-hoc estimation of BCI performance which is generally given in the BCI literature. For the fixed optimal condition, a re-sampling scheme (see section 2.5) is used to increase robustness.

Liu [17]. This approach simply sums the outputs of a linear classifier per class and over iterations. M is defined by

$\begin{equation} M(\boldsymbol{D},t,j) = \min _c\Bigg(\sum _{j^{\prime }=1}^j \boldsymbol{D}_{c,j^{\prime },t}\Bigg) \end{equation} \tag{ 4 }$

for trial t and iteration j, thus taking the most negative class. The average of all target classification scores in D^train, called $\tilde{d}$, is taken as the initialization for optimizing $\Delta = S \tilde{d}$ (section 2.5). Θ then contains this threshold Δ for online use.

During online operation, the decision for an early stop is taken as follows:

$\begin{equation} F(\boldsymbol{D}^{{\rm test}},\Theta ,t,j) = \left\lbrace \begin{array}{@{}ccc@{}}c_{j,t}^w & {\rm if}\, M(\boldsymbol{D}^{{\rm test}},t,j) < \Delta \\ \emptyset &\hspace{-41pt}{\rm otherwise} \end{array}\right. \end{equation} \tag{ 5 }$

where $c_{j,t}^w$ is the class with the most negative summed score.

In the original work, EEG features from iteration i are smoothed together by averaging them with features from iterations i − 1 and i − 2 before classification. Here, access is only given to the instantaneous classifier output by design, so this pre-processing step is left out.

Jin [20]. This method sets a minimum number of consecutive iterations (j^min) which must result in the same class prediction. To measure this, M(D, t, j) is defined as 1 + n, where n denotes the number of iterations directly prior to iteration j which all lead to the same decision as j does (according to equation (3)). The original paper performs an additional pre-processing step on D to give D* where $\boldsymbol{D}^*_{c,j,t} = \frac{1}{j} \sum ^{j}_{j^{\prime }=1}{\boldsymbol{D}_{c,j^{\prime },t}}$ for each j ⩾ j^mean. Since it only includes processing the classifier outputs, the same is done here. Both values of Θ = {j^min, j^mean} are themselves the hyperparameters (S₁ = j^min and S₂ = j^mean). They are optimized directly according to the procedure in section 2.5.

In the online case, the decision function F is as follows:

$\begin{equation} F(\boldsymbol{D}^{*,{\rm test}},\Theta ,t,j,) = \left\lbrace \begin{array}{@{}cc@{}}c_{j,t}^w & {\rm if}\, M(\boldsymbol{D}^{*,{\rm test}},t,j) \ge j^{{\rm min}} \\ & \hspace{-41pt}{\rm and}\, j \ge j^{{\rm mean}} \\ \emptyset , &\hspace{-52pt}{\rm otherwise.} \end{array}\right. \end{equation} \tag{ 6 }$

Here, both j^min and j^mean were bound between 2 and 5. As there are no decisions taken before j < j^mean, the earliest possible stop is at iteration j^mean + j^min.

Lenhardt [16]. This method was originally proposed for a visual speller, and defines two thresholds. For both thresholds, D is first transformed into D* by $\boldsymbol{D}^*_{1\ldots C,j,t} = {\rm norm}( \sum _{j^{\prime }=1}^j -\boldsymbol{D}_{1\ldots C,j^{\prime },t} )$, where ${\rm norm}(X) = \frac{X-\min (X)}{\max (X-\min (X))}$, a vector normalization over the class dimension. Lenhardt requires the target class to have positive labels. Since the target class was defined to be negative in this framework, a sign reversal is done within this transformation.

The first threshold is based on a measure called the 'intensity', or M₁, defined for iteration j and trial t as

$\begin{equation} M_{1}(\boldsymbol{D}^*,t,j) = \sum _{c=1}^{C} \boldsymbol{D}^*_{c,t,j} \end{equation} \tag{ 7 }$

M₁ is thus simply the sum of all values in D*_1...C,t,j. The second threshold is based on a slightly different measure M₂, defined as

$\begin{equation} M_2(\boldsymbol{D}^*,t,j) = 1-\frac{\boldsymbol{D}^{*}_{k_2,t,j}}{\boldsymbol{D}^{*}_{k_1,t,j}} \end{equation} \tag{ 8 }$

where $k_1 = {\rm argmax}_{c}(\boldsymbol{D}^*_{c,t,j})$ and $k_2 = {\rm argmax}_{c\ne k_1}(\boldsymbol{D}^*_{c,t,j})$. M₂ is thus the ratio of the highest to the second highest value in D*_1...C,t,j.

To find both thresholds, the training function L(D*^,train, l, Γ) first calculates the intensity M₁(D*, t, j) for the earliest possible stop which leads to a correct decision (according to equation (3)). This is done for each training trial, and the average of these intensities, $\tilde{i}$, serves as an initialization for optimizing $\Delta _1 = S_1 \tilde{i}$ (see section 2.5). The second threshold, Δ₂, is optimized on the training data directly (Δ₂ = S₂), which leads to Θ = {Δ₁, Δ₂}.

The online decisions for an early stop are taken by F as follows:

$\begin{eqnarray} \fl F(\boldsymbol{D}^{*,{\rm test}},\Theta ,t,j) = \left\lbrace \begin{array}{@{}cc@{}}c_{j,t}^w &{\rm if}\, M_1(\boldsymbol{D}^{*,{\rm test}},t,j) < \Delta _1 \\ & \hspace{7pt}{\rm and}\, M_2(\boldsymbol{D}^{*,{\rm test}},t,j) > \Delta _2 \\ \emptyset , & \hspace{-54pt}{\rm otherwise}. \end{array} \right. \end{eqnarray} \tag{ 9 }$

The original publication calculates the thresholds using data of several subjects. Here, it is restricted to the more common case of only having the calibration data of a single subject.

Rank diff [19]. Rank diff was designed for the AMUSE paradigm. However, the AMUSE dataset used in the current performance comparisons was not available during the design of Rank diff. The method is based on M(D, t, j), defined for iteration j and trial t as

$\begin{equation} M(\boldsymbol{D},t,j) = \min _{c^*}({\rm med}(\boldsymbol{D}_{c_{j,t}^w,1\ldots j,t})-{\rm med}(\boldsymbol{D}_{c^*,1\ldots j,t})) \end{equation} \tag{ 10 }$

where $c_{j,t}^w$ is found using equation (3), c* runs over all other classes, and med is the median operator. Θ = {Δ₁, ..., Δ_J} contains J iteration dependent thresholds Δ, allowing it to model the within-trial variation of ERP responses. To determine Δ_j (the threshold for iteration j), L(D^train, l, Γ, j) calculates M(D^train, t, j) for each training trial. These are assigned to one of two vectors, depending on whether the estimated class label was correct (+_j) or incorrect (−_j), using equation (3). To implement a security margin that prevents overly optimistic early stops, Δ_j is then defined as $\Delta _{j} = \max (\max ({\boldsymbol{-}}_j) \, , \, S\cdot {\rm med}({\boldsymbol{+}}_j))$. The hyperparameter S is set globally over all iterations, and is optimized as described in section 2.5.

A larger value for M means a better separation of the class medians, and thus more confidence in a decision. For online use, F is defined as

$\begin{equation} F(\boldsymbol{D}^{{\rm test}}, \Theta , t, j) = \left\lbrace \begin{array}{@{}cc@{}}c_{j,t}^w & {\rm if}\, M(\boldsymbol{D}^{{\rm test}},t,j) > \Delta _j \\ \emptyset ,&\hspace{-46pt}{\rm otherwise}. \end{array}\right. \end{equation} \tag{ 11 }$

Höhne [18]. This method was again proposed for an auditory speller paradigm. A t-test with unequal variances (also referred to as Welch's t-test) is calculated on D_1...j for each class against all other classes. Equation (12) shows the test statistic,

$\begin{equation} t(X_1, X_2) = \frac{\overline{X}_1 - \overline{X}_2}{\sqrt{\frac{s_1^2}{N_1} + \frac{s_2^2}{N_2}}} \end{equation} \tag{ 12 }$

where $\overline{X}_i, s_i^2, N_i$ are the sample mean, sample variance and sample size, respectively. In order to account for a varying number of data points, the t-test statistic $t(\boldsymbol{D}_{c,1\ldots j,t}, \boldsymbol{D}_{\tilde{c},1\ldots j,t})$ is transformed into a p-value (computed from a one-sided test, see supplementary appendix A (available from stacks.iop.org/JNE/10/036025/mmedia)). This transformation is part of any statistics toolbox. M(D, t, j) is then defined as this p-value.

The natural threshold on such a p-value is the significance level, or Θ = {α}. It is optimized for SPM on the training data directly, with α = S (section 2.5). During online operation, the decision for an early stop is taken according to the following

$\begin{equation} F(\boldsymbol{D}^{{\rm test}}, \Theta, t,j ) = \left\lbrace \begin{array}{@{}cc@{}}c_{j,t}^w &{\rm if}\, M(\boldsymbol{D}^{{\rm test}},t,j) < \alpha \\ \emptyset , &\hspace{-40pt}{\rm otherwise} \end{array}\right. \end{equation} \tag{ 13 }$

where $c_{j,t}^w$ is the class with the smallest M value.

Zhang [15]. This method proposes a probabilistic approach to continuously tracking whether a BCI user is in control or non-control state. The core of the algorithm is the formulation of a joint likelihood function, which accumulates the evidence for being attended for each stimulus over iterations. Here, this evidence is used as a criterion for early stopping.

The method requires that the distributions of the classifier outputs for target and non-target stimuli, f_a and f_u, are known. In practice, they are estimated from the training data. Here, univariate normal distributions for the classifier outputs are assumed, which are fully characterized by their means and variances μ_a/u and $\sigma ^2_{\rm a/u}$. Note that this assumption is compatible with the assumption of multivariate normal distributed data made by the LDA classifier (see supplementary appendix B (available from stacks.iop.org/JNE/10/036025/mmedia)).

Making the further simplifying assumption that the EEG data are independent across iterations, as well as across stimulations within one iteration, the joint likelihood of observing the classifier outputs of the first j iterations given that stimulus c is attended is

$\begin{equation} f\big({\boldsymbol{D}}^{{\rm test}}_{1 \dots C, 1 \dots j, t} \,\big|\, c\big) = \prod _{j^\ast =1}^{j} f_{\rm a}\big({\boldsymbol{D}}^{{\rm test}}_{c, j^\ast , t}\big) \prod _{c^\ast \ne c} f_{\rm u}\big({\boldsymbol{D}}^{{\rm test}}_{c^\ast , j^\ast , t}\big). \end{equation} \tag{ 14 }$

According to Bayes' theorem, the posterior probability of stimulus c being attended after the first j iterations is given by

$\begin{equation} P\big(c \,\big|\, {\boldsymbol{D}}^{{\rm test}}_{1 \dots C, 1 \dots j, t}\big) = \frac{f\big({\boldsymbol{D}}^{{\rm test}}_{1 \dots C, 1\dots j, t} \,\big|\, c\big) P(c)}{\sum _{\tilde{c} = 1}^C f\big({\boldsymbol{D}}^{{\rm test}}_{1 \dots c, 1 \dots j, t} \,\big|\, \tilde{c}\big) P(\tilde{c})} , \end{equation} \tag{ 15 }$

where P(c), c = 1, ..., C denote the prior probabilities with which the C classes are selected by the user. These probabilities may be estimated in the training phase or, e.g., derived from a language model in spelling applications. If all classes are assumed to be equally likely, uniform priors P(c) = 1/C, c = 1, ..., C are used. Note that in practice, numerical inaccuracies may occur when evaluating equations (14) and (15) due to excessively small values of the joint likelihood function. Therefore, all computations are carried out on log probabilities, as in the original publication [15].

The joint likelihoods and posterior probabilities of each class are updated after each iteration. The function M(D^test, t, j) is defined as the maximum of the posterior probabilities of all classes:

$\begin{equation} M(\boldsymbol{D}^{{\rm test}},t,j) = \max _{\check{c} = 1, \ldots , C} P\big(\check{c} \,\big|\, {\boldsymbol{D}}^{{\rm test}}_{1 \ldots C, 1 \ldots j, t}\big). \end{equation} \tag{ 16 }$

A natural stopping condition is attained, if equation (16) exceeds 1 − α, where α = S is a hyperparameter to be optimized based on the training data (section 2.5). Thus, $\Theta = \lbrace \lbrace P(c) \rbrace _{c = 1}^C, \mu _{\rm a/u}, \sigma ^2_{\rm a/u},\alpha \rbrace$. The decision at iteration j then reads:

$\begin{equation} F(\boldsymbol{D}^{{\rm test}}, \Theta ,t,j) = \left\lbrace \begin{array}{@{}cc@{}}c_{j,t}^w &{\rm if}\, M(\boldsymbol{D}^{{\rm test}},t,j) \ge 1 - \alpha \\ \emptyset ,&\hspace{-13pt}{\rm otherwise} \;\qquad \qquad , \end{array}\right. \end{equation} \tag{ 17 }$

where $c_{j,t}^w = {\rm argmax}_{\check{c} = 1, \ldots , C} P(\check{c} \,|\, {\boldsymbol{D}}^{{\rm test}}_{1 \ldots C, 1 \ldots j, t})$.

Upper bound. The upper bound method uses knowledge about the class labels and is not applicable for normal online use. It is applied for reasons of comparison and to deliver an upper bound on the performance of any of the other stopping methods. It simply stops a trial at the first occasion which leads to a correct decision. If a trial does not result in a correct decision at any iteration, the full number of iterations is applied.

The most important requirement for the methods is that their performance should at worst drop to the baseline performance of the no stopping strategy whenever data become difficult to separate. In other words, they should gracefully degrade. Several artificial datasets were constructed to test the robustness of the stopping methods for data distortions which are typically observed in classifier outputs during BCI online sessions (see figure 2). Unless stated otherwise, data are sampled from two normal distributions (target and non-target) with unit variance. They represent the classifier outputs and not the EEG data or features thereof. A training set and a validation set were generated independently. For the calculation of the SPM metric, all of the following datasets were assumed to have a SOA of 200 ms, six stimulus classes, a two-step interface, and an overhead of 10 s per step. The influence of the following distortion types were investigated (but not interactions between different types):

Zoom In Zoom Out Reset image size

Separation. To mimic data that is generally difficult to classify, the distance between the target and non-target class mean is varied from zero standard deviations (STD) for no separability, to three STD (good separability), effectively changing the AUC. This is done equally for both training and validation data. Based on the results, all other datasets were constructed with a relative high AUC value of 0.85.

Outliers. Artifacts in the EEG create outliers in the classifier outputs. These can generally be removed from the calibration data before training the classifier. Particularly during longer use, the number of artifacts in the online phase will be larger and removing them is not straightforward. Thus, between 0% (no outliers) and 65% (many outliers) are introduced to both classes, but in the validation data only. Outliers are randomly sampled from a uniform distribution on [−3, 3].

Confusions. A typical phenomenon in some BCIs is a systematical confusion, where one particular non-target is mistaken for a target more often than other non-targets [21]. This leads to a partial overlap of non-target and target classifiers scores. To simulate this, one non-target class was shifted from the non-target distribution (normal separation) toward the target distribution (no separation) in several steps. This is a special case of the separation dataset, where only one non-target stimulus class is affected by bad separability. It is applied to both training and validation data.

Datasets with changes in terms of drift, scale, distribution and number of classes were also tested. Since their influence on the performance was minimal, they are not described fully, but only depicted in figure 2.

EEG data from 83 sessions were taken, originating from AMUSE [19], an auditory ERP experimental BCI paradigm, and five visual ERP paradigms: rapid serial visual presentation (RSVP) [22], movement-onset visual evoked potential (MVEP) [23], and standard visual evoked potentials from the paradigms hex-o-spell, center speller and cake speller [10]. Datasets were previously recorded for the respective publications. All datasets consisted of an offline calibration part and an online spelling part; in all the original studies the users' task was to write text. Details on the datasets can be found in table 1. All calibration datasets were reduced to 24 trials. Due to the varying number of iterations per trial, there was still a discrepancy between datasets. Since it is not the primary goal here to compare datasets, this was not considered a problem. All online results are based on 21 training trials, to resemble the amount available in one cross-validation fold. Only results explicitly testing the influence of training data had a varying number of training trials.

2.3.1. Pre-processing and classification

For the evaluation and optimization of early stopping methods, a metric is needed which respects the speed–accuracy trade-off under consideration of the overhead. The information transfer rate (ITR) [25, 26] is an obvious performance measure, as it is predominant in the BCI literature. However, for spelling applications it is more straightforward and intuitive to use the SPM metric [27]. In addition to the overhead, it considers the application specific error handling, the benefit of which is aptly shown in figure 1. Typically, a correct selection counts as +1 symbol, and an erroneous selection counts as −1 symbol, to account for the added effort of a necessary subsequent corrective action. SPM is defined as follows:

$\begin{eqnarray} &&\fl {\mbox {effective correct}} = \mbox {percent correct - percent erroneous} \end{eqnarray} \tag{ 18 }$

$\begin{equation} {\rm SPM} = \mbox {decisions per minute} * {\mbox {effective correct}} \end{equation} \tag{ 19 }$

where decisions per minute considers the overhead and stimulation, as discussed in the introduction.

In a two-step selection process, errors can be made on the first level (group selection) and on the second level (single symbol selection). When the group selection is false, the second level contains the possibility to select an 'undo' action, which cancels the previous selection. As no symbol is selected after such a meander, it counts as a 0 symbol selection. Percent correct is now interpreted as the accuracy on both levels and percent erroneous as the percentage of erroneous selections on the second level. Note that the simplifying assumption that the accuracy is the same for each symbol was made (similarly to the assumption of the ITR as defined in [26]).

All methods have one or more hyperparameters S, which allow the user to influence the trade-off between speed and robustness. It requires some form of optimization within the training function L. Interestingly, most of the original publications do not provide a way to properly estimate the value of this hyperparameter.

Therefore, the hyperparameter S is optimized using a resampling procedure on the classifier scores obtained from a two-fold cross-validation on the training trials (see section 3). Two-thirds of the training trials are randomly selected to estimate the necessary data statistics (for instance the intensities in Lenhardt or the means, priors and variances in Zhang). Given these statistics, the remaining data is used to optimize hyperparameter S for SPM, using a simple linear search. This procedure is repeated 15 times and the median hyperparameter is subsequently applied to the test trials.

Some methods ((Fixed) optimal, Höhne, Jin) do not require such a compound optimization. However, they too are estimated using the resampling procedure for extra robustness and better comparability. In these cases the hyperparameter is optimized on the larger two-thirds of data split directly. In the case of two mutually dependent hyperparameters (Jin and Lenhardt), they are estimated simultaneously, by extending the linear search to a two dimensional grid search.

The size and heterogeneity of the datasets can be used to find proxies for the optimization of hyperparameter S. For this purpose, the optimized values of S for all training sets were investigated for dependencies on data features. Any resulting independent variable should be chosen such that it can easily be derived from a calibration dataset. As the AUC of the training data is highly correlated with S, and is independent of the sample size, it was chosen as the independent variable. The dependent variable S was then linearly fitted to the AUC of the training data. Fitting was done using iteratively re-weighted least-squares analysis [28] to reduce the influence of outliers. This gives two coefficients for each method, which transform the AUC of a dataset into a corresponding hyperparameter S†, which approximates S. All EEG datasets were re-analyzed using this approximation by S†, and the performance was compared to the optimized variant. In the case of multiple hyperparameters, the one with the highest correlation coefficient was approximated whilst fixing the second one. For Jin, j^mean was fixed to the minimum of two iterations, and for Lenhardt, S₁ was fixed to 1.

All of the methods are quite robust to most types of data distortion. The largest performance deteriorations were caused by the class separation, outliers and systematical confusions, which are shown in figure 4(B). Results are summarized by their mean performance gain over the no stopping condition (see figure 4(A)) over the given parameter space. The gain of a specific setting was considered significant if at least 90% of the 40 simulations yielded an output on the same side of the baseline (black bars). Bars are gray where no such significance was found.

Zoom In Zoom Out Reset image size

Separation. Most of the methods behave well for data that is difficult to separate, meaning that they do not decrease performance below the no stopping baseline when data separation goes to zero. Liu clearly struggles with the data as they becomes difficult to separate, as they severely undershoot the baseline performance. This is reflected in the relatively low average gain. To a lesser extent, Jin suffers from the same effect for the most difficult data. For difficult data, Rank diff shows a constant and tiny, but significant reduction. Zhang is exceptionally resilient to inseparable data, with a significant performance gain for all but the most inseparable data.

Outlier. All of the methods have severe problems in dealing with outliers. None of the methods reached a positive average gain over the given parameter space. Though Liu performs worst for outliers, here the picture is not that clear. In fact, all of the methods severely deteriorate performance compared to the baseline when many outliers are present, though the level of tolerance differs slightly between methods. It is important to note that the number of outliers went as high as 65%, and that these outliers were not observed in the calibration data.

Confusion. All methods are influenced by systematic confusions, but they degraded gracefully to the baseline. Interestingly, also Liu does not show any significant performance reduction from the baseline, although confusion is a special case of separation.

The results for all of the datasets can be found in the supplementary material (see supplementary figure 1 (available from stacks.iop.org/JNE/10/036025/mmedia)). It is worth noting that Liu is the only method that is clearly not robust against drifts, showing a performance peak around the no-drift parameter and falling off to both sides with positive and negative drifts. Furthermore, Liu, Rank diff, and Zhang show reduced performance gain when the data is downscaled. For all of the other methods, and for upscaling, only marginal influences were found. An increased number of classes reduced the gain for all methods, and a more skewed distribution showed higher gains for all methods.

While the results on the artificial data give an intuition of the properties of the various methods, the truly relevant question is their behavior on real BCI data, which is summarized in figure 5. On average over all datasets, all methods improve the performance. However, this does not hold when observing the paradigms individually.

Zoom In Zoom Out Reset image size

All of the methods significantly increase the performance on two datasets (center speller and RSVP speller), and all methods apart from fixed and fixed optimal managed to increase the performance on the cake speller dataset (for significance levels see figure 5). Performance on the hex-o-spell dataset is significantly increased by all but Liu. A significant increase in performance for the MVEP dataset is found for Rank diff, Zhang, and Lenhardt. None of the methods significantly increases (or decreases) performance on the AMUSE dataset.

In summary, no method increased performance significantly on all datasets. Rank diff, Zhang, and Lenhardt did so for five datasets, Höhne and Jin for four datasets, and fixed, fixed optimal, and Liu for three datasets.

Figure 6 shows the average gain over no stopping for each stopping method, accumulated over datasets. Over all datasets, Höhne, Lenhardt, and Zhang performed particularly well, as indicated by the high cumulative average gain (4.17, 4.14, and 4.47 respectively). Directly after these three methods, Jin, and Rank diff follow with cumulative gains of 3.11, and 3.49. The cumulative gain for Liu is competitive (2.83), in particular for the 'easier' datasets. When subtracting the negative gain for AMUSE and MVEP it reduces to 2.60.

Zoom In Zoom Out Reset image size

Both of the fixed methods (fixed and fixed optimal) show a cumulative gain over no stopping (2.08 and 2.36, after subtraction of negative paradigms). However, their gains are clearly smaller when compared to most dynamic methods. A full image depicting individual gains for all subjects and all stopping methods can be found in supplementary figure 2 (available from stacks.iop.org/JNE/10/036025/mmedia). In absolute terms, the fixed methods and Liu reduce performance for about 26.5% to 30% of the subjects, against 8.4% to 15.7% for the dynamic methods.

In short, Höhne, Lenhardt, and in particular Zhang perform best in terms of gain. Liu performs good for good subjects, but underperforms on difficult EEG data.

Figure 7 shows the effect of the training set size on the estimated online performance, as averaged over all datasets. With the minimum number of three training trials, all methods but Jin and Lenhardt decrease in performance with respect to no stopping, on average. With 12 training trials or more, all methods increase performance over the no stopping condition, on average. Performance for the no stopping condition converges after about 12 trials, as the accuracy saturates to 100% for several paradigms. For all stopping methods, the relative gain increases substantially during these first four conditions (three to 12 training trials). After 12 training trials, most dynamic methods keep increasing in performance until trial 21, showing a performance gain that is impossible without such stopping methods.

Zoom In Zoom Out Reset image size

Importantly, on average all dynamic methods, except for Liu, consistently outperform the fixed and fixed optimal method from six training trials upwards. Most dynamic approaches can thus exploit the data better than simply fixing a number of iterations, even with small amounts of training data. Liu needed on average 15 trials or more to outperform fixed (optimal). Though the effect size varies, the average findings largely transfer to the individual datasets (not shown).

The hyperparameter S of most methods is particularly correlated with the AUC of the training data. Even across different paradigms, with all the different settings, the optimized S is a rather robust function of this data feature (see supplementary figure 3 (available from stacks.iop.org/JNE/10/036025/mmedia)). The black line shows the fit of all of the data, and the gray lines fit the data, each with one dataset left out. Only for method Jin do these fits not clearly overlap, which is likely due to the fixing of j^mean in the online case, the inherent discrete nature of the hyperparameters and the small parameter range. All other methods have a robust fit, regardless of the aforementioned factors and the paradigm fitted. The resulting coefficients can be found in table 2 and can be used to approximate S according to the equation

$\begin{equation} S^\dagger = a\cdot {\rm AUC} + b. \end{equation} \tag{ 20 }$

The result of re-analyzing all of the data, using the above equation to determine S†, can be found in figure 8. Results are summarized by their median, 25 and 75 percentiles and minimum and maximum gain. Though the use of S† significantly decreased the training performance for several methods when compared to the optimal S (not shown), it did not significantly change the online performance for any method. Significance was tested using a Bonferroni-corrected, paired t-test.

Zoom In Zoom Out Reset image size

Early stopping methods allow the BCI to adapt to the current state of the user. In doing so, they can reduce the time required for a selection, and thus increase the robustness of the system. The latter seems to be particularly true for dynamic methods, since they can adapt to changes in the data. Though fixed and fixed optimal are the simplest methods and require very little adaptation of the BCI system, they are not dynamic; during online application the number of iterations is fixed regardless of the data coming in.

On average the fixed methods do increase the performance with respect to no stopping. Unfortunately, their performance gains leave much room for improvement, which is observed both in artificial data and real EEG data. This indeed is due to the fact that these methods are incapable of adapting to changes in the data distributions over time. When observing the online data, parts of the data that have an AUC below average show the largest reduction in performance for the fixed methods. Interestingly, they also show the largest performance improvement when non-stationarities temporarily lead to a higher SNR of the data.

As most of the original publications do not specify how to find an appropriate value for the hyperparameter (here called S), a general optimization scheme is proposed. With the SPM, it uses a realistic metric as objective function and leads to considerable performance gains. To optimize S, realistic simulations for the online performance for different values of S are required, which can be time consuming and error-prone to implement. For the BCI practitioner, a look-up table or regression function would be convenient. For instance, such a regression function would help to determine an at least suitable hyperparameter for a new BCI user, given some collected training data. The empirical observation that S is highly correlated with the AUC led to a set of coefficients that translate a given AUC of the training data into an approximation of S, called S†. Using S†, the online performance gains were similar to the performance gains using S.

Apart from an interesting basic finding, this result has practical relevance. The AUC can be estimated reliably with a relatively low amount of data, say a few trials. Using the regression coefficients, the problem with training on a low number of training trials may thus be overcome. But even with appropriate training size, they alleviate the practitioner of the task of running simulations to optimize S. They thereby reduce the barrier for using stopping methods in similar paradigms.

The classifier scores of the EEG data used in this study are available in the supplementary material (available from stacks.iop.org/JNE/10/036025/mmedia). It contains all data required for reproducing the current results and, more importantly, for benchmarking future methods. Classifier scores for each fold are given, such that over-fitting can be avoided and realistic comparisons can be made. Furthermore, example scripts were created to facilitate the use of these data. They take care of the proper cross-validation. Only one training and one test function need to be written for any new method.

Be aware that the implementation that is the basis for this paper differs slightly and contains a random factor; results will thus not exactly match those from the example script.

With offline analyses, the danger of overfitting the data is always tempting, and online validations are considered to be the last conclusive evidence. In this case, the protocol was kept strict so as to mimic such a validation, leaving the online data untouched until the very end. As running online trials with several stopping methods active at the same time would be impossible, a block structured protocol would be necessary. The current quasi-online protocol bypasses this need, whilst keeping the validity of the test intact. One clear difference between the two is that for the latter, the human is not in the loop. Potentially, shortening the trials can increase or reduce the workload on the user. In particular when an early stop results in a false decision, the user might get frustrated. This study shows the principle validity of the tested methods. Influences of such methods on the human in the loop remain an open topic, though at least for Jin [20], Lenhardt [16], and Rank diff [19] there is partial online evidence from direct comparison.

Method Liu differed slightly from the original implementation, to comply with the testing framework. The linear kernel support-vector machine [17] originally employed, is very similar to the LDA used here. This design difference is not likely to be the cause for the poor performance. The effect of temporally smoothing the features, as proposed by Liu, could result in cleaner data. At the same time it introduces strong dependencies between samples, which violates the assumptions of most classifiers. The influence of this on the current results is unclear.

For the sake of comparability, the current analyses were done using a single linear classifier. LDA is the most prominent classifier in the ERP based BCI literature, and particularly suited for ERP data [24]. However, given that for all methods the gained performance depends on the AUC of the classifier data, it could be expected that any improvement made to the AUC—by smarter feature extraction or better classification—could pay off double: it can improve the performance in itself, and additionally benefit more from the stopping methods. The influence of different, potentially nonlinear methods remains an open question.