Spatial Filtering for EEG-Based Regression Problems in Brain-Computer Interface (BCI)

TL;DR

This paper introduces two novel spatial filtering methods for EEG-based regression in brain-computer interfaces, significantly improving signal quality and estimation accuracy in a large-scale psychomotor vigilance study.

Contribution

The paper extends common spatial pattern (CSP) filters to regression problems using fuzzy sets, demonstrating their effectiveness in EEG response speed estimation.

Findings

Spatial filters significantly increase EEG signal quality.

Reduction of root mean square estimation error by up to 19.77%.

Increase in correlation to true response speed by up to 86.47%.

Abstract

Electroencephalogram (EEG) signals are frequently used in brain-computer interfaces (BCIs), but they are easily contaminated by artifacts and noises, so preprocessing must be done before they are fed into a machine learning algorithm for classification or regression. Spatial filters have been widely used to increase the signal-to-noise ratio of EEG for BCI classification problems, but their applications in BCI regression problems have been very limited. This paper proposes two common spatial pattern (CSP) filters for EEG-based regression problems in BCI, which are extended from the CSP filter for classification, by making use of fuzzy sets. Experimental results on EEG-based response speed estimation from a large-scale study, which collected 143 sessions of sustained-attention psychomotor vigilance task data from 17 subjects during a 5-month period, demonstrate that the two proposed spatial filters can significantly increase the EEG signal quality. When used in LASSO and k-nearest neighbors regression for user response speed estimation, the spatial filters can reduce the root mean square estimation error by 10.02-19.77%, and at the same time increase the correlation to the true response speed by 19.39-86.47%.

Tables

Table 1. TABLE I: p 𝑝 p -values of two-way ANOVA tests for { Raw , CAR , OVR , OVA } Raw CAR OVR OVA \{\texttt{Raw},\ \texttt{CAR},\ \texttt{OVR},\ \texttt{OVA}\} .

	LASSO		kNN
	RMSE	CC	RMSE	CC
$p$	$.0061$	$.0000$	$.0000$	$.0000$

Table 2. TABLE II: p 𝑝 p -values of non-parametric multiple comparison for { Raw , CAR , OVR , OVA } Raw CAR OVR OVA \{\texttt{Raw},\ \texttt{CAR},\ \texttt{OVR},\ \texttt{OVA}\} .

	LASSO						kNN
	RMSE			CC			RMSE			CC
	Raw	CAR	OVR	Raw	CAR	OVR	Raw	CAR	OVR	Raw	CAR	OVR
CAR	.5883			.3374			.1437			.0009
OVR	.0063	.0034		.0000	.0000		.0000	.0001		.0000	.0000
OVA	.0122	.0044	.4960	.0000	.0000	.4970	.0000	.0001	.4937	.0000	.0000	.4741

Equations

$$\displaystyle\mathbf{D}_{\tau_{k},1}=\left\{\begin{array}[]{ll}1,&\tau_{k}\mathrm{\ is\ the\ onset\ of\ the\ }\\ &k\mathrm{th\ target\ stimulus}\\ 0,&\mathrm{otherwise}\end{array}\right.$$

$${\mathbf{W}}^{\ast }=\arg \underset{\mathbf{W}}{\max }\frac{\mathrm{Tr}({\mathbf{W}}^T{\mathbf{PD}}^T{\mathbf{DP}}^T\mathbf{W})}{\mathrm{Tr}({\mathbf{W}}^T\mathbf{X}{\mathbf{X}}^T\mathbf{W})}$$

$${\mathbf{X}}_n^{\prime }={{\mathbf{W}}^{\ast }}^T{\mathbf{X}}_n,n=1,...,N.$$

$${\mathbf{W}}_k^{\ast }=\arg \underset{\mathbf{W}}{\max }\frac{\mathrm{Tr}({\mathbf{W}}^T{\bar{\mathbf{\Sigma }}}_k\mathbf{W})}{\mathrm{Tr}[{\mathbf{W}}^T({\sum }_{i\ne k}{\bar{\mathbf{\Sigma }}}_i)\mathbf{W}]}$$

$${\mathbf{W}}^{\ast }=[{\mathbf{W}}_1^{\ast },...{\mathbf{W}}_K^{\ast }]\in {\mathbb{R}}^{C\times KF}$$

$$A={\int }_{D_A}{\mu }_A(x)/x$$

$$p_k=\frac{100\cdot k}{K+1},k=1,...,K$$

$${\bar{\mathbf{X}}}_k=\frac{{\sum }_{n=1}^N{\mu }_k(y_n){\mathbf{X}}_n}{{\sum }_{n=1}^N{\mu }_k(y_n)},k=1,...,K$$

$${\mathbf{W}}_k^{\ast }=\arg \underset{\mathbf{W}}{\max }\frac{\mathrm{Tr}({\mathbf{W}}^T{\bar{\mathbf{\Sigma }}}_k\mathbf{W})}{\mathrm{Tr}[{\mathbf{W}}^T({\sum }_{i=1}^K{\bar{\mathbf{\Sigma }}}_i)\mathbf{W}]}$$

Full text

I Introduction

Electroencephalogram (EEG) signals are the most widely used input for brain-computer interfaces (BCIs) [24, 53, 25, 29, 47, 34], mainly due to the convenience to obtain them, compared with magnetoencephalography (MEG) [32], functional magnetic resonance imaging (fMRI) [44], functional near-infrared spectroscopy (fNIRS) [33], and invasive signals like electrocorticography (ECoG) [35] and intracortical neural recordings [30]. However, EEG signals are often contaminated by ocular, muscular, and cardiac artifacts and various noises (power-line, changes in electrode impedances, etc) [49, 4, 34]. Usually some preprocessing, either manually or automatically [4, 34], is needed to remove the artifacts, and then temporal and spatial filters are applied to further improve the EEG signal quality before feeding it into a classification or regression algorithm. The most commonly used temporal filters are band-pass filters and notch filters (at 50 or 60 Hz power-line frequency).

This paper focuses on spatial filtering for improving the EEG signal quality. Many such approaches have been proposed in the literature [54, 17, 38, 40, 41, 15, 7, 37, 2]. However, almost all of them focus primarily on EEG classification problems in BCI, whereas EEG regression problems have been largely overlooked. Nevertheless, the latter is also very important in BCI. One example is driver drowsiness (or alertness) estimation from EEG signals, which has been extensively studied in our previous research [59, 27, 26, 62, 60, 28, 56]. This is a very important problem because drowsy driving is among the most important causes of road crashes, following only to alcohol, speeding, and inattention [43]. According to the National Highway Traffic Safety Administration [52], 2.5% of fatal motor vehicle crashes (on average 886/year in the U.S.) and 2.5% of fatalities (on average 1,004/year in the U.S.) between 2005 and 2009 involved drowsy driving.

This paper proposes two spatial filters for EEG-based regression problems in BCI. We also validate their performance in response speed (RS) estimation from EEG signals measured in a large-scale sustained-attention psychomotor vigilance task (PVT) [21], which collected 143 sessions of data from 17 subjects in a 5-month period.

The remainder of this paper is organized as follows: Section II reviews the state-of-the-art spatial filters for EEG-based classification problems in BCI. Section III introduces our proposed spatial filters for supervised BCI regression problems. Section IV describes the experimental setup, RS and EEG data preprocessing techniques, and the procedure to evaluate the performances of different spatial filters. Section V presents the results of the comparative studies and parameter sensitivity analysis for the proposed spatial filter. Section VI discusses the limitations of the proposed approaches and outlines several future research directions. Finally, Section VII draws conclusions.

II Spatial Filters for EEG Classification in BCI

Many spatial filters have been proposed for EEG classification in BCI. The most basic ones include common average reference (CAR) [48], Laplacian filters [23], and principal component analysis [19]. Some of the more recent and also more sophisticated ones are:

Independent Component Analysis (ICA) [9, 54, 17], which decomposes a multivariate signal into independent non-Gaussian signals. ICA has been widely used in the EEG research community to detect and remove stereotyped eye, muscle, and line noise artifacts [20, 49, 26].

Generally ICA works on an unepoched long block of EEG data, instead of epoched short EEG trials. Let the unepoched EEG data be $\mathbf{X}\in\mathbb{R}^{C\times T}$, where $C$ is the number of EEG channels, and $T$ is the number of time samples. ICA assumes that $\mathbf{X}$ is the linear combination of $c$ independent sources, i.e., $\mathbf{X}=\mathbf{AS}$, where $\mathbf{A}\in\mathbb{R}^{C\times c}$ is the mixing matrix, and the source signals, which are the rows of $\mathbf{S}\in\mathbb{R}^{c\times T}$, are supposed to be stationary, independent, and non-Gaussian. ICA can use various different principles [49, 17, 54, 9] to estimate both unknown $\mathbf{A}$ and unknown $\mathbf{S}$ simultaneously from $\mathbf{X}$. Once $\mathbf{S}$ is obtained, cleaner and more representative features may be extracted from it than from the original $\mathbf{X}$ [26]. 2. 2.

xDAWN algorithm [38, 39, 40], which is often used to increase the signal to signal-plus-noise ratio in P300-based BCIs.

Like ICA, xDAWN also works on the unepoched long block of EEG data $\mathbf{X}\in\mathbb{R}^{C\times T}$. It assumes that $\mathbf{X}=\mathbf{PD}^{T}+\mathbf{N}$, where $\mathbf{P}\in\mathbb{R}^{C\times S}$ represents the P300 signal in an EEG epoch, and $\mathbf{D}\in\mathbb{R}^{T\times S}$ is a Toeplitz matrix whose first column is defined as:

[TABLE]

and $\mathbf{N}\in\mathbb{R}^{C\times T}$ represents the ongoing background brain activity as well as the artifacts and noises. xDAWN then designs a spatial filtering matrix $\mathbf{W}^{*}\in\mathbb{R}^{C\times F}$, where $F$ is the number of spatial filters, to maximize the signal to signal-plus-noise ratio, i.e.,

[TABLE]

where $\mathrm{Tr}(\cdot)$ is the trace of a matrix. (5) is a generalized Rayleigh quotient [14], and its solution $\mathbf{W}^{*}$ is the concatenation of the $F$ eigenvectors associated with the $F$ largest eigenvalues of the matrix $(\mathbf{X}\mathbf{X}^{T})^{-1}\mathbf{PD}^{T}\mathbf{DP}^{T}$.

The spatially filtered trial for $\mathbf{X}_{n}$ is then computed as:

[TABLE] 3. 3.

Canonical Correlation Analysis (CCA) [41, 15], which finds linear transformations to maximize the correlations between two datasets. It has been used to improve BCI performance in code-modulated visual evoked potentials [5], steady-state visual evoked potentials [6], and event-related potentials like P300 and error-related potentials [45].

Unlike ICA and xDAWN, CCA works on epoched EEG trials. Consider a binary classification problem, with $N_{1}$ training examples in Class 1 and $N_{2}$ training examples in Class 2. Let $(\mathbf{X}_{n},y_{n})$ be the $n$th training example, where $\mathbf{X}_{n}\in\mathbb{R}^{C\times S}$ ($C$ is the number of channels, and $S$ is the number of time samples in each trial), and $y_{n}\in\{1,2\}$. Let $\bar{\mathbf{X}}_{k}\in\mathbb{R}^{C\times S}$ be the average of $\mathbf{X}_{n}$ in Class $k$ $(k=1,2)$. We then construct $\tilde{\mathbf{X}}=[\tilde{\mathbf{X}}_{1}\ \tilde{\mathbf{X}}_{2}]$ and $\tilde{\mathbf{Z}}=[\tilde{\mathbf{Z}}_{1}\ \tilde{\mathbf{Z}}_{2}]$, where $\tilde{\mathbf{X}}_{k}$ is the concatenation of all $N_{k}$ $\mathbf{X}_{n}$ in Class $k$, and $\tilde{\mathbf{Z}}_{k}$ is the concatenation of $N_{k}$ $\bar{\mathbf{X}}_{k}$. CCA first finds two vector filters $\mathbf{w}_{\tilde{\mathbf{X}}}$ and $\mathbf{w}_{\tilde{\mathbf{Z}}}$ such that the correlation between $\mathbf{w}_{\tilde{\mathbf{X}}}^{T}\tilde{\mathbf{X}}$ and $\mathbf{w}_{\tilde{\mathbf{Z}}}^{T}\tilde{\mathbf{Z}}$ is maximized. $\mathbf{w}_{\tilde{\mathbf{X}}}^{T}\mathbf{X}$ and $\mathbf{w}_{\tilde{\mathbf{Z}}}^{T}\tilde{\mathbf{Z}}$ are called the first pair of canonical variables. CCA then finds the second pair of canonical variables in a similar way, subject to the constraint that they are uncorrelated with the first pair of canonical variables. This procedure can be continued up to $C$ times.

Finally, the spatial filtering matrix is the concatenation of all $\mathbf{w}_{\tilde{\mathbf{X}}}$, which can be applied to each $\mathbf{X}_{n}$ to increase its SNR. 4. 4.

Common Spatial Patterns (CSP) [7, 37], which is a supervised technique frequently used to enhance the binary classification performance of EEG data. The basic idea is to separate the EEG signal into additive subcomponents which have maximum differences in variance between the two classes. In the following we introduce the one-versus-the-rest (OVR) CSP [11], which extends the traditional CSP from binary classification to $K$ classes.

Like CCA, OVR CSP also works on epoched EEG trials. Let $(\mathbf{X}_{n},y_{n})$ be the $n$th training example, as defined above. Assume the mean of $\mathbf{X}_{n}$ has been removed, e.g., by high-pass or band-pass filtering. Then, for Class $k$, OVR CSP finds a spatial filter matrix $\mathbf{W}_{k}^{*}\in\mathbb{R}^{C\times F}$, where $F$ is the number of spatial filters, to maximize the variance difference between Class $k$ and the rest:

[TABLE]

where $\bar{\mathbf{\Sigma}}_{k}$ is the mean covariance matrix of trials in Class $k$. (7) is also a generalized Rayleigh quotient [14], and the solution $\mathbf{W}_{k}^{*}$ is the concatenation of the $F$ eigenvectors associated with the $F$ largest eigenvalues of the matrix $(\sum_{i\neq k}\bar{\mathbf{\Sigma}}_{i})^{-1}\bar{\mathbf{\Sigma}}_{k}$.

Finally, we concatenate the $K$ individual OVR CSP spatial filters to obtain the complete filter:

[TABLE]

and compute the spatially filtered trial for $\mathbf{X}_{n}$ by (6).

III Spatial Filters for Supervised BCI Regression Problems

In this section we propose two common spatial pattern for regression (CSPR) filters, which extend the multi-class CSP filters from classification to regression by making use of fuzzy sets [63], as we have done in [62].

First, a brief introduction of fuzzy sets is given below.

III-A Fuzzy Sets

A fuzzy set $A$ is comprised of a universe of discourse $D_{A}$ of real numbers together with a membership function $\mu_{A}:\ D_{A}\to[0,1]$, i.e.,

[TABLE]

Here $\int$ denotes the collection of all points $x\in D_{A}$ with associated membership degree $\mu_{A}(x)$. An example of a fuzzy set is shown in Fig. 1. The membership degrees are $\mu_{A}(1)=0$, $\mu_{A}(3)=0.5$, $\mu_{A}(5)=1$, $\mu_{A}(6)=0.8$, and $\mu_{A}(10)=0$. Observe that this is different from traditional (binary) sets, where each element can only belong to a set completely (i.e., with membership degree 1), or does not belong to it at all (i.e., with membership degree 0); there is nothing in between (i.e., with membership degree 0.5). Fuzzy sets are frequently used in modeling concepts in natural language [22, 55, 36], which may not have clear boundaries.

III-B CSPR-OVR

Let $\mathbf{X}_{n}\in\mathbb{R}^{C\times S}$ ($n=1,...,N$) be the $n$th EEG trial, where $C$ is the number of channels and $S$ is the number of time samples in each trial. We assume that the mean of each channel measurement has been removed, which is usually performed by band-pass filtering. Let $y_{n}\in\{1,...,K\}$ be the RS of $\mathbf{X}_{n}$.

With the help of fuzzy sets, we can define “fuzzy” classes to connect regression problems and classification problems. Assume $K$ fuzzy classes are used. First, we partition the interval $[0,100]$ into $K+1$ equal intervals, and denote the partition points as $\{p_{k}\}_{k=1,...,K}$. It is easy to obtain that

[TABLE]

For each $p_{k}$, we then find the corresponding $p_{k}$ percentile value of all training $y_{n}$ and denote it as $P_{k}$. Next we define $K$ fuzzy classes from them, as shown in Fig. 2. In this way, we can “classify” the training $y_{n}$ into $K$ fuzzy classes, corresponding to the $K$ crisp classes in the CSP for classification. However, note that in the CSP for classification a $y_{n}$ belongs to a crisp class either completely or not at all. For a fuzzy class here, a $y_{n}$ can belong to it at a membership degree in $[0,1]$.

Next, for each fuzzy class, we compute its mean EEG trial as:

[TABLE]

where $\mu_{k}(y_{n})$ is the membership degree of $y_{n}$ in Fuzzy Class $k$. Substituting (11) into (7), we can solve for the spatial filtering matrix $\mathbf{W}_{k}^{*}$ for Fuzzy Class $k$. Essentially, this $\mathbf{W}_{k}^{*}$ makes those $\mathbf{X}_{n}$ in Fuzzy Class $k$ different from those not in Fuzzy Class $k$, which will help the regression performance, as we will demonstrate in Section V.

Next, we construct the concatenated spatial filtering matrix $\mathbf{W}^{*}$ by (8), and finally perform the spatial filtering for each EEG trial $\mathbf{X}_{n}$ by (6). The complete CSPR-OVR spatial filter for supervised BCI regression problems is summarized in Algorithm 1.

III-C CSPR-OVA

In (7) we construct the multi-class CSP using an OVR approach, but it can also be constructed using the following one-versus-all (OVA) approach:

[TABLE]

The only difference between (12) and (7) is that the numerator of (12) also includes the contribution from Class $k$ itself. If we view Class $k$ as the signal of interest, and all other classes as noises, then (12) maximizes the signal to signal-plus-noise ratio, as (5) in the xDAWN algorithm.

Equation (12) is also a generalized Rayleigh quotient [14], and the solution $\mathbf{W}_{k}^{*}$ is the concatenation of the $F$ eigenvectors associated with the $F$ largest eigenvalues of the matrix $(\sum_{i=1}^{K}\bar{\mathbf{\Sigma}}_{i})^{-1}\bar{\mathbf{\Sigma}}_{k}$. The OVA CSP for classification still uses (8) to construct the final spatial filter, and (6) to perform the filtering.

Using the technique introduced in the previous subsection, we can easily develop the CSPR-OVA spatial filter for BCI regression problems. Its procedure is almost identical to that in Algorithm 1. The only difference is that $\mathbf{W}_{k}^{*}$ is computed by (12) instead of (7).

IV Experiments and Data

This section introduces a PVT experiment that was used to evaluate the performances of the proposed spatial filtering algorithms, the corresponding RS and EEG data preprocessing procedures, and the feature sets.

IV-A Experiment Setup

17 university students (13 males; average age 22.4, standard deviation 1.6) from National Chiao Tung University (NCTU) in Taiwan volunteered to support the data-collection efforts over a 5-month period to study EEG correlates of attention and performance changes under specific conditions of real-world fatigue [21], as determined by the effectiveness score of Readiband [42]. The voluntary, fully informed consent of the persons used in this research was obtained as required by federal and Army regulations [51, 50]. The Institutional Review Board of NCTU approved the experimental protocol.

All participants registered their fatigue levels through a smartphone daily, and received notifications to report for experimental trials when the effectiveness score deemed their conditions fitted the experimental requirement (low fatigue: $>90$; normal: $[70,90]$; high fatigue: $<70$). Upon completion of the related questionnaires [Karolinska Sleepiness Scale (KSS) [1], and electronically-adapted visual analog scale for fatigue (VAS-F) and stress (VAS-S)] and the informed consent form, subjects performed a PVT, a dynamic attention-shifting task, a lane-keeping task, and selected surveys (KSS, VAS-F, VAS-S, state-trait anxiety inventory, and mind-wandering) preceding each condition. EEG data were recorded at 1000 Hz using a 64-channel NeuroScan system. Most participants performed the laboratory experiment thrice in each of the three fatigue states.

In this paper we focus on the PVT [10], which is a sustained-attention task that uses RS to measure the speed with which a subject responds to a visual stimulus. It is widely used, particularly by NASA, for its ease of scoring, simple metrics, convergent validity, and free of learning effects. In our experiment, the PVT was presented on a smartphone with each trial initiated as an empty solid white circle centered on the touchscreen that began to fill in red displayed as a clockwise sweeping motion like the hand of a clock. The sweeping motion was programmed to turn solid red in one second or terminate upon a response by the participants, which required them to tap the touchscreen with the thumb of their dominant hand. The RS was computed as the inverse of the elapsed time between the appearance of the empty solid white circle and the participant’s response. Following completion of each trial, the circle went back to solid white until the next trial. Inter-trial intervals consisted of random intervals between 2-10 seconds.

143 sessions of PVT data were collected from the 17 subjects, and each session lasted 10 minutes. Our goal is to predict the RS using a 3-second EEG trial immediately before it.

IV-B Performance Evaluation Process

The following procedure was performed to evaluate the performances of different spatial filters:

EEG data preprocessing to suppress artifacts and noises. 2. 2.

RS data preprocessing to suppress outliers. 3. 3.

5-fold cross-validation to compute the regression performance for each combination of spatial filters and regression method: first randomly partition the trials into five equal folds; then, use four folds for supervised spatial filtering and regression model training, and the remaining fold for testing; repeat this five times so that every fold is used in testing; finally compute the regression performances in terms of root mean square error (RMSE) and correlation coefficient (CC). Two regression methods were used: LASSO, whose adjustable parameter $\lambda$ was optimized by an inner 5-fold cross-validation on the training dataset, and $k$-nearest neighbors (kNN) regression, where $k=5$. 4. 4.

Repeat Step 3 10 times and compute the average regression performance.

More details about the first two steps are given in the next two subsections.

IV-C EEG Data Preprocessing

We first downsampled the EEG data to 256 Hz, then epoched them to 3-second trials according to the onset of the PVTs. Let the onset time of the $n$th PVT be $t_{n}$. Then, the 62-channel EEG trial in $[t_{n}-3,t_{n}]$ seconds was used to predict the RS, i.e., $\mathbf{X}_{n}\in\mathbb{R}^{62\times 768}$. Each trial was then individually filtered by a $[1,20]$ Hz finite impulse response band-pass filter to make each channel zero-mean and to remove un-useful high frequency components.

Because the inter-trial intervals consisted of random intervals between 2-10 seconds, it’s possible that a 3-second EEG trial covers part of data from the previous trial. Additionally, a trial may also contain the EEG oscillations related to motor reaction (tapping the touchscreen) in the previous trial. To remedy these problems, we removed overlapping trials: let the RS of the $n$th trial be $y_{n}$ (the corresponding response time is $1/y_{n}$); then, the $n$th trial is removed if $t_{n}-t_{n-1}<1/y_{n-1}+3$, i.e., when the 3-second EEG data for Trial $n$ overlap with the data and response for the previous trial.

IV-D RS Data Preprocessing

The raw response times for two subjects are shown in Fig. 3. The top panel is from a typical subject, whose response times were mostly shorter than 1 second. The lower panel is from a subject with possible data recording issues, because lots of response times were longer than 5 seconds, which are highly unlikely in practice. So we excluded that subject from consideration in this paper, and only used the remaining 16 subjects.

As shown in Fig. 3, the response times were very noisy, and there were obvious outliers. It is very important to suppress the outliers and noises so that the performances of different algorithms can be more accurately compared. In addition to the step in the previous subsection to remove overlapping trials, we also employed the following 2-step procedure for response time preprocessing:

Outlier thresholding, which aimed to suppress abnormally large response times. First, a threshold $\theta=m_{y}+3\sigma_{y}$ was computed for each subject, where $m_{y}$ is the mean response time from all sessions of that subject, and $\sigma_{y}$ is the corresponding standard deviation. Then, all response times larger than $\theta$ were replaced by $\theta$. Note that the threshold was different for different subjects. 2. 2.

Moving average smoothing, which replaced each response time by the average response time during a 60 seconds moving window centered at the onset of the corresponding PVT to suppress noises.

We then computed the RS as the inverse of the RT. The RSs for the 16 subjects are shown in Fig. 4. Observe that they are roughly in the same range, and many of them are approximately Gaussian.

IV-E Feature Extraction

We extracted the following four feature sets for each preprocessed EEG trial:

• •

Raw: Theta and Alpha powerband features from the band-pass filtered EEG trials. We computed the average power spectral density (PSD) in the Theta band (4-8 Hz) and Alpha band (8-13 Hz) for each channel using Welch’s method [57], and converted these $62\times 2=124$ band powers to dBs as our features.

• •

CAR: Theta and Alpha powerband features from EEG trials filtered by CAR. This procedure was almost identical to Raw, except that the band-pass filtered EEG trials were also spatially filtered by CAR before the $62\times 2=124$ powerband features were computed. CAR is one of the most commonly used spatial filters for EEG, and [31] showed that it helped improve EEG classification performance. It simply removes the mean of all channels from each channel.

• •

OVR: Theta and Alpha powerband features from EEG trials filtered by CSPR-OVR. This procedure was almost identical to CAR, except that the CAR filter was replaced by CSPR-OVR. We used 3 fuzzy classes for the RSs, and 21 spatial filters111We used 21 spatial filters here so that the filtered signals had roughly the same dimensionality as the original signals, which ensured fair performance comparison. In Section V-C we also performed sensitivity analysis on the number of spatial filters. for each fuzzy class, so that the spatially filtered signals had dimensionality $63\times 1280$, roughly the same as the dimensionality of the original signals. We then extracted $63\times 2=126$ band power features for each trial.

• •

OVA: Theta and Alpha powerband features from EEG trials filtered by CSPR-OVA. This procedure was also almost identical to CAR, except that the spatial filtering was performed by CSPR-OVA instead of CAR. There were also $63\times 2=126$ band power features for each trial.

V Experimental Results

This section compares the informativeness of the features in Raw, CAR, OVR and OVA, presents the regression performances, and also performs parameter sensitivity analysis for Algorithm 1.

V-A Informativeness of the Features

Before studying the regression performances, it is important to check if the extracted features in Raw, CAR, OVR and OVA are indeed meaningful. We picked a typical subject, partitioned his data random into 50% training and 50% testing, and extracted Raw and CAR. We then designed the spatial filters using CSPR-OVR and CSPR-OVA on the training data, and extracted the corresponding OVR and OVA. For each feature set, we identified the top three channels that had the maximum correlation with the RS using the training data, and also computed the corresponding correlation coefficients for the testing data.

The results are shown in Fig. 5, where in each subfigure the data on the left of the black dotted line were used for training, and the right for testing. The top thick curve is the RS, and the bottom three curves are the maximally correlated features (note that good features are negatively correlated with the RS) identified from the training data. The training and testing correlation coefficients are shown on the left and right of the corresponding channel, respectively. Observe that the features from CAR had slightly better correlations with the RS in training than those from Raw, but not necessarily in testing. However, the features from OVR and OVA had much higher training and testing correlations to the RS than those from Raw and CAR, suggesting that CSPR-OVR and CSPR-OVA can indeed increase the signal quality. The reason is: if we view Class $k$ as the signal of interest, and all other classes as noises, then CSPR-OVR in (7) enhances the signal to noise ratio of the EEG signal, and CSPR-OVA in (12) enhances the signal to signal-plus-noise ratio.

V-B Regression Performance Comparison

The RMSEs and CCs of LASSO and kNN using the four feature sets are shown in Fig. 6 for the 16 subjects. Recall that for each subject the feature extraction methods were run 10 times, each with randomly partitioned training and testing data, and the average regression performances are shown here. The average RMSEs and CCs across all subjects are also shown in the last group of each panel. Observe that CAR had comparable or slightly better performance than Raw. Regardless of which regression algorithm was used, generally OVR and OVA had similar performance, and both of them achieved much smaller RMSEs and much larger CCs than Raw and CAR, suggesting that our extension of CSP from supervised classification to supervised regression can indeed improve the regression performance. Finally, LASSO had better performance than kNN on Raw and CAR, but kNN became better on OVR and OVA.

The corresponding percentage performance improvements of LASSO and kNN using the four feature sets are shown in Fig. 7, where the legend “LASSO,OVR/Raw” means the percentage performance improvement of LASSO on OVR over LASSO on Raw, and other legends should be interpreted in a similar manner. For both LASSO and kNN, OVR and OVA achieved similar performance improvements over Raw, and also over CAR. For LASSO, on average OVR had $10.02\%$ smaller RMSE than Raw, and $19.39\%$ larger CC. For kNN, on average OVR had $19.77\%$ smaller RMSE than Raw, and $86.47\%$ larger CC.

We also performed a two-way Analysis of Variance (ANOVA) for different regression algorithms to check if the RMSE and CC differences among the four feature sets were statistically significant, by setting the subjects as a random effect. The results are shown in Table I, which indicated that there were statistically significant differences in both RMSEs and CCs among different feature sets for both LASSO and kNN.

Then, non-parametric multiple comparison tests based on Dunn’s procedure [12, 13] were used to determine if the difference between any pair of algorithms was statistically significant, with a $p$-value correction using the False Discovery Rate method [3]. The $p$-values are shown in Table II, where the statistically significant ones are marked in bold. Table II shows that, except for the CC of kNN, generally there was no statistically significant difference between Raw and CAR. However, for both LASSO and kNN, the RMSE and CC differences between $\{\texttt{OVR},\ \texttt{OVA}\}$ and $\{\texttt{Raw},\ \texttt{CAR}\}$ were always statistically significant. In all cases, there were no statistically significant differences between OVR and OVA.

V-C Parameter Sensitivity Analysis

There are two adjustable parameters in CSPR-OVR: $K$, the number of fuzzy classes for the RSs, and $F$, the number of spatial filters for each fuzzy class. In this subsection we study the sensitivity of the regression performance to these two parameters.

The regression performances for $K=\{2,3,4,5,6,7\}$ ($F$ was fixed to be 21) are shown in Fig. 8. Algorithm 1 was repeated five times, each with a random partition of training and testing data, and the average regression results are shown. For both LASSO and kNN, on average $K=2$ gave worst performance, but $K=\{3,4,5,6,7\}$ resulted in roughly the same RMSE and CC. Hence, $K=3$ seems to be a good compromise between performance and computational cost.

The regression performances for $F=\{5,10,20,30,40,50,60\}$ ($K$ was fixed to be 3) are shown in Fig. 9. Algorithm 1 was again repeated five times, and the average regression results are shown. For both LASSO and kNN, generally a larger $F$ resulted in a smaller RMSE and a larger CC, but the performance may reach a plateau at a certain $F$. Also, a larger $F$ means heavier computational cost, which should be taken into consideration in choosing $F$. For the PVT experiment, $F\in[20,30]$ seemed to achieve a good compromise between performance and computational cost.

V-D Different Fuzzy Set Shapes

In Section III we used triangular fuzzy sets for simplicity, but other shapes can also be used. Fig. 10 illustrates how Gaussian fuzzy sets can be designed here: the center of the $k$th Gaussian fuzzy class is at $P_{k}$ [computed from (10)], and the spread is specially designed so that two adjacent fuzzy sets intersect at the midpoint with membership grade 0.5. As a result, generally the Gaussian fuzzy classes are not symmetric.

When the Gaussian fuzzy classes in Fig. 10 are used in CSPR-OVR and CSPR-OVA, the results are shown in Fig. 11, which are almost identical to those obtained from triangular fuzzy sets (Fig. 6).

V-E Robustness to Noise

It is also important to study the robustness of different spatial filters to noises. According to [64], there are two types of noises: class noise, which is the noise on the model outputs, and attribute noise, which is the noise on the model inputs. In this subsection we focus on the attribute noise.

As in [64], for each model input, we randomly replaced $q\%$ ($q=0,\,10,...,40$) of all trials from a subject with a uniform noise between its minimum and maximum values. After this was done for both the training and testing data, we extracted feature sets Raw, CAR, OVR and OVA, and trained LASSO and kNN, on the corrupted training data. We then tested their performances on the corrupted testing data. The results are shown in Fig. 12. Generally, as the noise level increased, the performances decreased, which is intuitive. However, OVR and OVA achieved better RMSEs and CCs than Raw and CAR at almost all noise levels, suggesting that it is still beneficial to use CSPR-OVR and CSPR-OVA even under high attribute noise.

V-F Computational cost

Observe from Algorithm 1 that in training CSPR-OVR needs to perform a matrix inversion and an eigen-decomposition to compute $\mathbf{W}^{*}$; however, once the training is done, the filtering of new EEG trials can be conducted very efficiently by a simple matrix multiplication [see (6)]. Let $N$ be the number of training samples. Then, the actual training time of CSPR-OVR and CSPR-OVA increased linearly with $N$, as shown in Fig. 13. The platform was a Dell XPS15 laptop (Intel i7-6700HQ CPU @2.60GHz, 16 GB memory) running Windows 10 Pro 64-bit and Matlab 2016b. A least squares curve fit shows that the training time is $0.2216+0.0003N$ seconds, which should not be a problem for a practical $N$.

VI Discussions and Future Research

Recall that 5-fold cross-validation was used in the performance evaluation in the previous section, i.e., we concatenated the nine-session data from the same subject, randomly partitioned them into five equal-length folds, and then used four folds for training and the remaining one for testing. So, the training and testing folds contained data from the same sessions. This is equivalent to the case that we label some session-specific data in offline regression. Our results showed that in this case CSPR-OVR and CSPR-OVA can significantly improve the regression performance.

To avoid the use of session-specific data, we also investigated a different validation method: leave-one-session-out validation, in which for each subject we trained the spatial filters using eight sessions and tested them on the remaining session. Interestingly, all four feature sets and both regression models achieved very poor performance here. The reasons are: 1) we need a proper way to normalize the RSs from different sessions, as done for the response times in [16]; and, 2) there is large intra-subject variation, meaning that the EEG responses for the same subject vary at different times (recall that these nine sessions were collected at different days); so, the patterns learned from previous sessions become obsolete for the new session, and hence spatial filtering alone does not help. However, our previous research [58, 62, 61] has shown that transfer learning can cope well with the inter-subject variation (individual differences) in both classification and regression problems, and we conjecture that it can also handle the intra-subject variation. One of our future research directions is to demonstrate the performance of CSPR-OVR and CSPR-OVA in a transfer learning framework to individualize a generalized model for regression problems, as done in [18, 46] for EEG-based cognitive performance classification.

Another direction of our future research will apply CSPR-OVR and CSPR-OVA to other important EEG-based regression problems, e.g., drowsiness (or alertness) estimation during driving, and integrate it with more sophisticated feature extraction approaches, e.g., Riemannian geometry [8], for better regression performance.

VII Conclusions

EEG signals are easily contaminated by artifacts and noises, so preprocessing is needed before they are fed into a machine learning algorithm in BCI. Spatial filters, e.g., ICA, xDAWN, CSP and CCA, have been widely used to increase the EEG signal quality for classification problems, but their applications in BCI regression problems have been very limited. In this paper, we have proposed two CSP filters for EEG-based regression problems in BCI, which were extended from the CSP filter for classification, by making use of fuzzy sets. Extensive experimental results on EEG-based RS estimation from a large-scale study, which collected 143 sessions of PVT data from 17 subjects during a 5-month period, demonstrated that our proposed spatial filters can significantly increase the EEG signal quality. When used in LASSO and kNN, the spatial filters can reduce the estimation RMSE by $10.02-19.77\%$, and at the same time increase the CC by $19.39-86.47\%$.

Acknowledgement

Research was sponsored by the U.S. Army Research Laboratory and was accomplished under Cooperative Agreement Numbers W911NF-10-2-0022 and W911NF-10-D-0002/TO 0023. The views and the conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Army Research Laboratory or the U.S. Government. This work was also partially supported by the Australian Research Council (ARC) under discovery grant DP150101645.