Driver Drowsiness Estimation from EEG Signals Using Online Weighted Adaptation Regularization for Regression (OwARR)

TL;DR

This paper introduces OwARR, a novel online transfer learning algorithm that improves driver drowsiness estimation from EEG signals with minimal calibration data, demonstrating superior accuracy and efficiency.

Contribution

It presents OwARR, a new online weighted adaptation regularization method for regression in BCI, integrating fuzzy sets and source domain selection to reduce calibration effort and computational cost.

Findings

OwARR achieves significantly lower estimation errors than existing methods.

OwARR-SDS reduces computational cost by about half.

The methods demonstrate robustness across subjects.

Abstract

One big challenge that hinders the transition of brain-computer interfaces (BCIs) from laboratory settings to real-life applications is the availability of high-performance and robust learning algorithms that can effectively handle individual differences, i.e., algorithms that can be applied to a new subject with zero or very little subject-specific calibration data. Transfer learning and domain adaptation have been extensively used for this purpose. However, most previous works focused on classification problems. This paper considers an important regression problem in BCI, namely, online driver drowsiness estimation from EEG signals. By integrating fuzzy sets with domain adaptation, we propose a novel online weighted adaptation regularization for regression (OwARR) algorithm to reduce the amount of subject-specific calibration data, and also a source domain selection (SDS) approach to save about half of the computational cost of OwARR. Using a simulated driving dataset with 15 subjects, we show that OwARR and OwARR-SDS can achieve significantly smaller estimation errors than several other approaches. We also provide comprehensive analyses on the robustness of OwARR and OwARR-SDS.

Tables

Table 1. TABLE I: p 𝑝 p -values of non-parametric multiple comparisons.

	OwARR	OwARR	OwARR	OwARR-SDS	OwARR-SDS	OwARR-SDS	OwARR-SDS
$m$	vs BL1	vs BL2	vs DAMF	vs BL1	vs BL2	vs DAMF	vs OwARR
0	.0007	N/A	.4073	.0011	N/A	.5091	.5000
5	.0000	.0000	.0000	.0000	.0000	.0000	.4813
10	.0000	.0000	.0000	.0000	.0000	.0000	.3888
15	.0000	.0000	.0000	.0000	.0000	.0000	.4075
20	.0000	.0000	.0000	.0000	.0000	.0000	.4795
25	.0000	.0000	.0000	.0000	.0000	.0000	.4658
30	.0000	.0000	.0000	.0000	.0000	.0000	.4153
35	.0000	.0000	.0001	.0000	.0000	.0001	.4364
40	.0000	.0000	.0002	.0000	.0000	.0004	.3956
45	.0000	.0000	.0004	.0000	.0000	.0006	.4378
50	.0000	.0000	.0005	.0000	.0000	.0013	.3766
55	.0000	.0000	.0007	.0000	.0000	.0035	.2898
60	.0000	.0000	.0011	.0000	.0000	.0108	.2132
65	.0000	.0000	.0015	.0000	.0000	.0097	.2573
70	.0000	.0000	.0024	.0000	.0000	.0143	.2527
75	.0000	.0000	.0038	.0000	.0000	.0133	.3122
80	.0000	.0000	.0047	.0000	.0000	.0294	.2304
85	.0000	.0000	.0058	.0000	.0000	.0513	.1785
90	.0000	.0000	.0071	.0000	.0000	.0383	.2387
95	.0000	.0000	.0086	.0000	.0000	.1091	.1205
100	.0000	.0000	.0117	.0000	.0000	.1179	.1352

Equations

$$f(\mathbf{x})=Q(y|\mathbf{x})=\frac{P(\mathbf{x},y)}{P(\mathbf{x})}=\frac{Q(\mathbf{x}|y)P(y)}{P(\mathbf{x})},$$

$$f(\mathbf{x})={\mathbf{\alpha }}^T\mathbf{x}$$

$$f=$$

$$+\lambda [d(P^z,P^t)+d(Q^z,Q^t)]-\gamma {\tilde{r}}^2(y,f(\mathbf{x}))$$

$$X$$

$$\mathbf{y}$$

$$\displaystyle E_{ii}=\left\{\begin{array}[]{ll}1,&1\leq i\leq n\\ w^{t},&n+1\leq i\leq n+m\end{array}\right.$$

$$\sum _{i=1}^n(y_i-f({\mathbf{x}}_i){)}^2+w^t\sum _{i=n+1}^{n+m}(y_i-f({\mathbf{x}}_i){)}^2$$

$$=$$

$$w^t=\max (2,\sigma \cdot n/m)$$

$$d(P^z,P^t)$$

$$={\mathbf{\alpha }}^{T}X{M}_{P}X\mathbf{\alpha }$$

$$\displaystyle(M_{P})_{ij}=\left\{\begin{array}[]{ll}\frac{1}{n^{2}},&1\leq i\leq n,1\leq j\leq n\\ \frac{1}{m^{2}},&n+1\leq i\leq n+m,\\ &n+1\leq j\leq n+m\\ \frac{-1}{nm},&\text{otherwise}\end{array}\right.$$

$$d(Q^z,Q^t)=\sum _{c=1}^C{\left[\frac{1}{n_c}\sum _{{\mathbf{x}}_i\in {\mathcal{D}}_c^z}f({\mathbf{x}}_i)-\frac{1}{m_c}\sum _{{\mathbf{x}}_i\in {\mathcal{D}}_c^t}f({\mathbf{x}}_i)\right]}^2$$

$${\bar{\mu }}_{ic}^z=\frac{{\mu }_{ic}^z}{{\sum }_{i=1}^n{\mu }_{ic}^z},i=1,...,n;c=1,2,3$$

$$d(Q^z,Q^t)=\sum _{c=1}^3{\left[\sum _{{\mathbf{x}}_i\in {\mathcal{D}}^z}{\bar{\mu }}_{ic}^{z}f({\mathbf{x}}_i)-\sum _{{\text{x}}_i\in {\mathcal{D}}^t}{\bar{\mu }}_{ic}^{t}f({\mathbf{x}}_i)\right]}^2$$

$$d(Q^z,Q^t)$$

$$=\sum _{c=1}^3{\mathbf{\alpha }}^{T}X{M}_{c}X\mathbf{\alpha }={\mathbf{\alpha }}^{T}X{M}_{Q}X\mathbf{\alpha }$$

$$M_Q=M_1+M_2+M_3$$

$$\displaystyle(M_{c})_{ij}=\left\{\begin{array}[]{ll}\bar{\mu}_{ic}^{z}\bar{\mu}_{jc}^{z},&\mathbf{x}_{i},\,\mathbf{x}_{j}\in\mathcal{D}_{c}^{z}\\ \bar{\mu}_{ic}^{t}\bar{\mu}_{jc}^{t},&\mathbf{x}_{i},\,\mathbf{x}_{j}\in\mathcal{D}_{c}^{t}\\ -\bar{\mu}_{ic}^{z}\bar{\mu}_{jc}^{t},&\mathbf{x}_{i}\in\mathcal{D}_{c}^{z},\,\mathbf{x}_{j}\in\mathcal{D}_{c}^{t}\\ -\bar{\mu}_{ic}^{t}\bar{\mu}_{jc}^{z},&\mathbf{x}_{i}\in\mathcal{D}_{c}^{t},\,\mathbf{x}_{j}\in\mathcal{D}_{c}^{z}\\ 0,&\mbox{otherwise}\end{array}\right.$$

$$r(y,f(\mathbf{x}))$$

$$=\frac{{\mathbf{y}}^{T}X\mathbf{\alpha }}{\sqrt{{\mathbf{y}}^T\mathbf{y}}\cdot \sqrt{{\mathbf{\alpha }}^T{X}^{T}X\mathbf{\alpha }}}$$

$$r^2(y,f(\mathbf{x}))=\frac{{\mathbf{\alpha }}^T{X}^T{\mathbf{yy}}^{T}X\mathbf{\alpha }}{{\mathbf{y}}^T\mathbf{y}\cdot {\mathbf{\alpha }}^T{X}^{T}X\mathbf{\alpha }}$$

$${\tilde{r}}^2(y,f(\mathbf{x}))$$

$$=\frac{{\mathbf{\alpha }}^T{X}^T({\mathbf{yy}}^T-I)X\mathbf{\alpha }}{{\mathbf{y}}^T\mathbf{y}}$$

$$\mathbf{\alpha }=$$

$$+\lambda {\mathbf{\alpha }}^T{X}^T(M_P+M_Q)X\mathbf{\alpha }$$

$$+\gamma \frac{{\mathbf{\alpha }}^T{X}^T(I-{\mathbf{yy}}^T)X\mathbf{\alpha }}{{\mathbf{y}}^T\mathbf{y}}$$

$$\mathbf{\alpha }={\left[X^T\left(E+\lambda M_P+\lambda M_Q+\gamma \frac{I-{\mathbf{yy}}^T}{{\mathbf{y}}^T\mathbf{y}}\right)X\right]}^{-1}{X}^{T}E\mathbf{y}$$

$$d(z,t)=\sum _{c=1}^3||{\mathbf{m}}_c^z-{\mathbf{m}}_c^t||$$

$$y=\max \left\{0,\frac{1-e^{-(\tau -{\tau }_0)}}{1+e^{-(\tau -{\tau }_0)}}\right\}$$

$$X={\int }_{D_X}{\mu }_{{}_X}(x)/x$$

Full text

I Introduction

Brain computer interfaces (BCIs) [32, 45, 18, 28, 53] have attracted rapidly increasing research interest in the last decade, thanks to recent advances in neurosciences, wearable/mobile biosensors, and analytics. However, there are still many challenges in their transition from laboratory settings to real-life applications, including the reliability and convenience of the sensing hardware [21], and the availability of high-performance and robust algorithms for signal analysis and interpretation that can effectively handle individual differences and non-stationarity [25, 28, 50, 12]. This paper focuses on the last challenge, more specifically, how to generalize a BCI algorithm to a new subject, with zero or very little subject-specific calibration data.

Transfer learning (TL) [34], which improves learning in a new task by leveraging data or knowledge from other relevant tasks, represents a promising solution to the above challenge. Many TL approaches have been proposed for BCI applications [50], including: 1) feature representation transfer [13, 7, 39, 41], which encodes the knowledge across different tasks as features; 2) instance transfer [19, 20, 63, 56], which uses certain parts of the data from other tasks to help the learning for the current task; and, 3) classifier transfer, which includes domain adaptation (DA) [47, 1, 41], ensemble learning [43, 42], and their combinations [61, 60, 54].

However, most of the above TL approaches consider only BCI classification problems. Reducing the calibration data requirement in BCI regression problems has been largely under-studied. One example is online driver drowsiness estimation from EEG signals, which will be investigated in this paper. 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 [38]. According to the National Highway Traffic Safety Administration [44], 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. However, to our best knowledge, there have been only two works [51, 54] on TL for drowsiness estimation. Wei et al. [51] showed that selective TL, which selectively turns TL on or off based the level of session generalizability, can achieve better estimation performance than approaches that always turn TL on or off. Wu et al. [54] proposed a domain adaptation with model fusion (DAMF) approach for drowsiness estimation. By making use of data from other subjects in a DA framework, DAMF requires very little subject-specific calibration data, which significantly increases its real-world applicability.

In this paper, by making use of fuzzy sets (FSs) [67], we extend our earlier work on online weighted adaptation regularization [60] from classification to regression to estimate driver drowsiness online from EEG signals. We show that our two proposed algorithms can achieve significantly better estimation performance than the DAMF and two other baseline approaches.

The remainder of the paper is organized as follows: Section II introduces the details of the proposed online weighted adaptation regularization for regression (OwARR) algorithm. Section III further introduces a source domain selection (SDS) approach to save the computational cost of OwARR. Section IV presents experimental results and performance comparisons of OwARR and OwARR-SDS with three other approaches. Finally, Section V draws conclusions and points out future research directions.

II Online Weighted Adaptation Regularization for Regression (OwARR)

In [60] we have defined two types of calibration in BCI:

Offline calibration, in which a pool of unlabeled EEG epochs have been obtained a priori, and a subject or an oracle is queried to label some of these epochs, which are then used to train a model to label the remaining epochs in the pool. 2. 2.

Online calibration, in which some labeled EEG epochs are obtained on-the-fly, and then a model is trained from them for future (unseen) EEG epochs.

The major different between them is that, for offline calibration, the unlabeled EEG epochs can be used to help design the model (e.g., semi-supervised learning), whereas in online calibration there are no unlabeled EEG epochs. Additionally, in offline calibration we can query any epoch in the pool for the label, but in online calibration usually the sequence of the epochs is pre-determined and the subject or oracle has little control on which epochs to see next.

We only consider online calibration in this paper. This section introduces the OwARR algorithm, which extends the online weighted adaptation regularization algorithm [60] from classification to regression, by making use of FSs.

II-A Problem Definition

A domain [34, 23] $\mathcal{D}$ in TL consists of a $d$-dimensional feature space $\mathcal{X}$ and a marginal probability distribution $P(\mathbf{x})$, i.e., $\mathcal{D}=\{\mathcal{X},P(\mathbf{x})\}$, where $\mathbf{x}\in\mathcal{X}$. Two domains $\mathcal{D}^{z}$ and $\mathcal{D}^{t}$ are different if $\mathcal{X}^{z}\neq\mathcal{X}^{t}$, and/or $P^{z}(\mathbf{x})\neq P^{t}(\mathbf{x})$.

A task [34, 23] $\mathcal{T}$ in TL consists of an output space $\mathcal{Y}$ and a conditional probability distribution $Q(y|\mathbf{x})$. Two tasks $\mathcal{T}^{z}$ and $\mathcal{T}^{t}$ are different if $\mathcal{Y}^{z}\neq\mathcal{Y}^{t}$, or $Q^{z}(y|\mathbf{x})\neq Q^{t}(y|\mathbf{x})$.

Given the $z^{\mathrm{th}}$ source domain $\mathcal{D}^{z}$ with $n_{z}$ samples $(\mathbf{x}_{i}^{z},y_{i}^{z})$, $i=1,...,n_{z}$, and a target domain $\mathcal{D}^{t}$ with $m$ calibration samples $(\mathbf{x}_{j}^{t},y_{j}^{t})$, $j=1,...,m$, DA aims to learn a target prediction function $f(\mathbf{x}):\mathbf{x}\mapsto y$ with low expected error on $\mathcal{D}^{t}$, under the assumptions that $\mathcal{X}^{z}=\mathcal{X}^{t}$, $\mathcal{Y}^{z}=\mathcal{Y}^{t}$, $P^{z}(\mathbf{x})\neq P^{t}(\mathbf{x})$, and $Q^{z}(y|\mathbf{x})\neq Q^{t}(y|\mathbf{x})$.

In driver drowsiness estimation from EEG signals, EEG signals from a new subject are in the target domain, while EEG signals from the $z^{\mathrm{th}}$ existing subject are in the $z^{\mathrm{th}}$ source domain. A single data sample consists of the feature vector for a single EEG epoch in either domain. Though the features in source and target domains are extracted in the same way, generally their marginal and conditional probability distributions are different, i.e., $P^{z}(\mathbf{x})\neq P^{t}(\mathbf{x})$ and $Q^{z}(y|\mathbf{x})\neq Q^{t}(y|\mathbf{x})$, because different subjects usually have similar but distinct drowsy neural responses. As a result, data from a source domain cannot represent data in the target domain accurately, and must be integrated with some target domain data to induce the target domain regression function.

II-B The Learning Framework

Because

[TABLE]

to use the data in the $z^{\mathrm{th}}$ source domain in the target domain, we need to minimize the distance between the marginal and conditional probability distributions in the two domains by ensuring that111Strictly speaking, we should also make sure $P^{z}(y)$ is close to $P^{t}(y)$. In this paper we assume $P^{z}(y)$ and $P^{t}(y)$ are close. Our future research will consider the general case that $P^{z}(y)$ and $P^{t}(y)$ are different. $P^{z}(\mathbf{x})$ is close to $P^{t}(\mathbf{x})$, and $Q^{z}(\mathbf{x}|y)$ is also close to $Q^{t}(\mathbf{x}|y)$.

Assume both the output and each dimension of the input vector have zero mean. Then, the regression function can be written as

[TABLE]

where $\boldsymbol{\alpha}$ is the regression parameter vector to be found. The learning framework of OwARR is then formulated as:

[TABLE]

where $\lambda$ and $\gamma$ are non-negative regularization parameters, and $w^{t}$ is the overall weight for target domain samples, which should be larger than 1 so that more emphasis is given to target domain samples than source domain samples.

Briefly speaking, the first two terms in (3) minimize the sum of squared errors in the source domain and target domain, respectively. The 3rd term minimizes the distance between the marginal and conditional probability distributions in the two domains. The last term maximizes the approximate sample Pearson correlation coefficient between $y$ and $f(\mathbf{x})$, which helps avoid the undesirable situation that the regression output is (nearly) a constant.

In the next subsections we will explain how to compute the individual terms in (3).

II-C Sum of Squared Error Minimization

Let

[TABLE]

where the first $n$ $\mathbf{x}_{i}$ and $y_{i}$ are the column input vectors and the corresponding outputs in the source domain, the next $m$ $\mathbf{x}_{i}$ and $y_{i}$ are the column input vectors and the corresponding outputs in the target domain, and $T$ is the matrix transpose operation.

Define $E\in R^{(n+m)\times(n+m)}$ as a diagonal matrix with

[TABLE]

Then, the first two terms in (3) can be rewritten as

[TABLE]

The optimal selection of $w^{t}$ is very important to the performance of the OwARR algorithm. In this paper we use

[TABLE]

where $\sigma$ is a positive adjustable parameter, based on the following heuristics: 1) when $m$ is small, each target domain sample should have a large weight so that the target domain is not overwhelmed by the source domain; 2) as $m$ increases, the weight on the target domain samples should decrease gradually so that the source domain is not overwhelmed by the target domain; and, 3) the target domain samples should always have larger weights than the source domain samples because eventually the regression model will be applied to the target domain.

II-D Marginal Probability Distribution Adaptation

As in [23, 36, 60, 61], we compute $d(P^{z},P^{t})$ using the maximum mean discrepancy (MMD):

[TABLE]

where $M_{P}\in R^{(n+m)\times(n+m)}$ is the MMD matrix:

[TABLE]

II-E Conditional Probability Distribution Adaptation

In [23, 60, 61] a classification problem is considered, and it is more straightforward to perform conditional probability distribution adaptation. In this subsection we first briefly introduce the technique used there, and then describe in detail how we can perform conditional probability distribution adaptation in regression in a similar way, with the help of FSs [67, 16], which have been widely used in EEG feature extraction [5, 24, 15] and pattern recognition [33, 22, 11].

II-E1 Conditional Probability Distribution Adaptation for Classification

Let $\mathcal{D}_{c}^{z}=\{\mathbf{x}_{i}|\mathbf{x}_{i}\in\mathcal{D}^{z}\wedge y_{i}=c\}$ be the set of samples in Class $c$ ($c=1,...,C$) of the $z^{\mathrm{th}}$ source domain, $\mathcal{D}_{c}^{t}=\{\mathbf{x}_{i}|\mathbf{x}_{i}\in\mathcal{D}^{t}\wedge y_{i}=c\}$ be the set of samples in Class $c$ of the target domain, $n_{c}=|\mathcal{D}_{c}^{z}|$, and $m_{c}=|\mathcal{D}_{c}^{t}|$. Then, the distance between the conditional probability distributions in source and target domains is computed as the sum of the Euclidian distances between the class means in the two domains [23, 61, 60], i.e.,

[TABLE]

II-E2 Conditional Probability Distribution Adaptation for Regression

With the help of FSs (background materials are given in Appendix), we can transform the regression problem into a “classification” problem and hence perform conditional probability distribution adaptation using (17). First, for the $n_{z}$ outputs, $\{y_{i}^{z}\}_{i=1,...,n}$, in the $z^{\mathrm{th}}$ source domain, we find their 5, 50 and 95 percentile222There is a popular regression analysis method called quantile regression [17] in statistics and econometrics, which estimates either the conditional median or other quantiles of the response variable. The percentiles used in this paper are found directly from the data, and they should not be confused with quantile regression. values, $p_{5}^{z}$, $p_{50}^{z}$ and $p_{95}^{z}$, respectively, and define three triangular FSs333There can be other ways to define these FSs, e.g., we could use Gaussian FSs instead of triangular FSs, use $p_{10}^{z}$, $p_{50}^{z}$ and $p_{90}^{z}$ instead of $p_{5}^{z}$, $p_{50}^{z}$ and $p_{95}^{z}$, use other than three FSs in each domain, or use type-2 FSs [30] instead of type-1. Three type-1 triangular FSs are used here for simplicity. More discussions on the sensitivity of the OwARR algorithm to the number of type-1 triangular FSs are given in Section IV-H., $Small^{z}$, $Medium^{z}$ and $Large^{z}$, based on them, as shown in Fig. 1. In this way, we can “classify” the outputs in the $z^{\mathrm{th}}$ source domain into three fuzzy classes, $Small^{z}$, $Medium^{z}$ and $Large^{z}$, corresponding to the different classes in a traditional crisp classification problem. However, note that in the traditional crisp classification problem a sample can only belong to one class. For the fuzzy classes here, a sample can belong to more than one class simultaneously, at different degrees.

Denote Class $Small^{z}$ as Class 1, Class $Medium^{z}$ as Class 2, Class $Large^{z}$ as Class 3, and the membership degree of $y_{i}^{z}$ in Class $c$ as $\mu_{ic}^{z}$. We then normalize each $\mu_{ic}^{z}$ according to its class, i.e.,

[TABLE]

Similarly, we also find $p_{5}^{t}$, $p_{50}^{t}$ and $p_{95}^{t}$ from the $m$ target domain outputs $\{y_{i}^{t}\}_{i=n+1,...,n+m}$, define three FSs, $Small^{t}$, $Medium^{t}$ and $Large^{t}$, as shown in Fig. 1, and compute the corresponding normalized $\bar{\mu}_{ic}^{t}$, $i=n+1,...,n+m$, $c=1,2,3$.

Finally, similar to (17), the distance between the conditional probability distributions in the target domain and the $z^{\mathrm{th}}$ source domain is computed as:

[TABLE]

Substituting (2) into (19), it follows that

[TABLE]

where

[TABLE]

in which $M_{1}$, $M_{2}$ and $M_{3}$ are MMD matrices computed as:

[TABLE]

II-F Maximize the Approximate Sample Pearson Correlation Coefficient

The sample Pearson correlation coefficient $r(y,f(\mathbf{x}))$ is defined as [48]:

[TABLE]

and hence

[TABLE]

Note that $r^{2}(y,f(\mathbf{x}))$ has $\boldsymbol{\alpha}$ in the denominator, so it is very challenging to find a closed-form solution to maximize it. However, observe that $r^{2}(y,f(\mathbf{x}))$ increases as $\boldsymbol{\alpha}^{T}X^{T}\mathbf{yy}^{T}X\boldsymbol{\alpha}$ increases, and decreases as $\boldsymbol{\alpha}^{T}X^{T}X\boldsymbol{\alpha}$ increases. So, instead of maximizing $r^{2}(y,f(\mathbf{x}))$ directly, in this paper we try to maximize the following function:

[TABLE]

where $I\in R^{(n+m)\times(n+m)}$ is an identity matrix. $\tilde{r}^{2}(y,f(\mathbf{x}))$ has the same property as $r^{2}(y,f(\mathbf{x}))$, i.e., $\tilde{r}^{2}(y,f(\mathbf{x}))$ increases as $\boldsymbol{\alpha}^{T}X^{T}\mathbf{yy}^{T}X\boldsymbol{\alpha}$ increases, and decreases as $\boldsymbol{\alpha}^{T}X^{T}X\boldsymbol{\alpha}$ increases.

II-G The Closed-Form Solution

Substituting (9), (11), (20) and (30) into (3), we can rewrite it as

[TABLE]

Setting the derivative of the objective function above to 0 leads to

[TABLE]

II-H The Complete OwARR Algorithm

The pseudo-code for the complete OwARR algorithm is described in Algorithm 1. We first perform OwARR for each source domain separately, and then construct the final regression model as a weighted average of these base models, where the weight is the inverse of the training accuracy of the corresponding base model. The final regression model will then be applied to future unlabeled data.

III OwARR-SDS

A SDS procedure for online classification problems has been proposed in [60]. In this paper it is extended to regression problems by using the fuzzy classes again defined in Fig. 1. The primary goal of SDS is to reduce the computational cost of OwARR, because when there is a large number of source domains, performing OwARR for each source domain and then aggregating the base models would be very time-consuming.

Assume there are $Z$ different source domains. For the $z^{\mathrm{th}}$ source domain, we first compute $\mathbf{m}_{c}^{z}$ ($c=1,2,3$), the mean vector of each fuzzy class. Then, we also compute $\mathbf{m}_{c}^{t}$, the mean vector of each fuzzy class in the target domain, from the $m$ labeled samples. The distance between the two domains is:

[TABLE]

We next cluster the $Z$ numbers, $\{d(z,t)\}_{z=1,...,Z}$, by $k$-means clustering, and finally choose the cluster that has the smallest centroid, i.e., the source domains that are closest to the target domain. In this way, on average we only need to perform OwARR for $Z/k$ source domains. We used $k=2$ in this paper.

The pseudo-code for the complete OwARR-SDS algorithm is described in Algorithm 2. We first use SDS to select the $Z^{\prime}$ closest source domains, and then perform DA for each selected source domain separately. The final regression model is a weighted average of these base models, with the weight being the inverse of the training accuracy of the corresponding base model.

IV Experiments and Discussions

Experimental results on driver drowsiness estimation from EEG signals are presented in this section to demonstrate the performance of OwARR and OwARR-SDS.

IV-A Experiment Setup

We reused the experiment setup and data in [54]. 16 healthy subjects with normal or corrected to normal vision were recruited to participate in a sustained-attention driving experiment [4, 3], which consisted of a real vehicle mounted on a motion platform with 6 degrees of freedom immersed in a 360-degree virtual-reality scene. Each participant read and signed an informed consent form before the experiment began. Each experiment lasted for about 60-90 minutes and was conducted in the afternoon when the circadian rhythm of sleepiness reached its peak. To induce drowsiness during driving, the virtual-reality scenes simulated monotonous driving at a fixed 100 km/h speed on a straight and empty highway. During the experiment, lane-departure events were randomly applied every 5-10 seconds, and participants were instructed to steer the vehicle to compensate for these perturbations as quickly as possible. Subjects’ cognitive states and driving performance were monitored via a surveillance video camera and the vehicle trajectory throughout the experiment. The response time in response to the perturbation was recorded and later converted to drowsiness index. Meanwhile, participants’ scalp EEG signals were recorded using a 32-channel (30-channel EEGs plus 2-channel earlobes) 500 Hz Neuroscan NuAmps Express system (Compumedics Ltd., VIC, Australia).

The Institutional Review Board of the Taipei Veterans General Hospital approved the experimental protocol.

IV-B Evaluation Process and Performance Measures

The complete procedure for the application of OwARR for driver drowsiness estimation is shown in Algorithm 3. Compared with Algorithm 1 on OwARR for a generic application, here we also include detailed EEG data pre-processing and feature extraction steps. Observe that although the feature extraction methods for different auxiliary subjects have the same steps, their parameters (channels removed, principal components used, ranges used in normalization) may be different, so we need to record them for each auxiliary subject so that the features in the individual regression models can be computed correctly.

From the experiments we already knew the drowsiness indices for all $\sim$1200 epochs. To evaluate the performances of different algorithms, for each subject, we used up to 100 epochs in a randomly chosen continuous block for calibration, and the rest $\sim$1100 epochs for testing. Every time when five epochs were acquired, we computed the testing performance to show how the performance of the regression models changed over time. We ran this evaluation process 30 times, each time with a randomly chosen 100-epoch calibration block, to obtain statistically meaningful results. Finally, we repeated this entire process 15 times so that each subject had a chance to be the “15th” subject.

The primary performance measured used in this paper is the root mean squared error (RMSE) between the $\sim$1100 true drowsiness indices and the corresponding estimates for the testing epochs, which is optimized in the object functions of all algorithms. The secondary performance measure is the correlation coefficient (CC) between the true drowsiness indices and the estimates.

IV-C Preprocessing and Feature Extraction

The 16 subjects had different lengths of experiment, because the disturbances were presented randomly every 5-10 seconds. Data from one subject was not correctly recorded, so we used only 15 subjects. To ensure fair comparison, we used only the first 3,600 seconds data for each subject.

We defined a function [51, 54] to map the response time $\tau$ to a drowsiness index $y\in[0,1]$:

[TABLE]

$\tau_{0}=1$ was used in this paper, as in [54]. The drowsiness indices were then smoothed using a 90-second square moving-average window to reduce variations. This does not reduce the sensitivity of the drowsiness index because the cycle lengths of drowsiness fluctuations are longer than 4 minutes [26]. The smoothed drowsiness indices for the 15 subjects are shown in Fig. 2. Observe that each subject had some drowsiness indices at or close to 1, indicating drowsy driving.

We used EEGLAB [6] for EEG signal preprocessing. A band-pass filter (1-50 Hz) was applied to remove high-frequency muscle artifacts, line-noise contamination and DC drift. Next the EEG data were downsampled from 500 Hz to 250 Hz and re-referenced to averaged earlobes.

We tried to predict the drowsiness index for each subject every three seconds. All 30 EEG channels were used in feature extraction. We epoched 30-second EEG signals right before each sample point, and computed the average power spectral density (PSD) in the theta band (4-7.5 Hz) for each channel using Welch’s method [52], as research [27] has shown that theta band spectrum is a strong indicator of drowsiness. The theta band powers for three selected channels and the corresponding drowsiness index for Subject 1 are shown in Fig. 3. Observe that drowsiness index has strong correlations with the theta band powers.

Next, we converted the 30 theta band powers to dBs. To remove noises or bad channel readings, we removed channels whose maximum dBs were larger than 20. We then normalized the dBs of each remaining channel to mean zero and standard deviation one, and extracted a few (usually around 10) leading principal components, which accounted for 95% of the variance. The projections of the theta band powers onto these principal components were then normalized to $[0,1]$ and used as our features. Three such features for Subject 1 are shown in Fig. 3. Observe that the score on the first principal component has obvious correlation with the drowsiness index, suggesting that estimating the drowsiness index from the scores on the principal components is possible.

IV-D Algorithms

We compared the performances of OwARR and OwARR-SDS with three other algorithms introduced in [54]:

Baseline 1 (BL1), which combines data from all 14 existing subjects, builds a ridge regression model [10], and applies it to the new subject. That is, BL1 tries to build a subject-independent regression model and ignores data from the new subject completely. 2. 2.

Baseline 2 (BL2), which builds a ridge regression model using only subject-specific calibration samples from the new subject. That is, BL2 ignores data from existing subjects completely. 3. 3.

DAMF, which builds 14 ridge regression models by combining data from each auxiliary subject with data from the new subject, respectively, and then uses a weighted average to obtain the final regression model. The weights are also the inverse of the training RMSEs, as in Algorithms 1 and 2.

The ridge parameter $\sigma=0.01$ was used in the above three algorithms, as in [54]. For OwARR and OwARR-SDS, we used $\sigma=0.2$, $\lambda=10$, and $\gamma=0.5$. However, as will be shown in Section IV-H, OwARR and OwARR-SDS are robust to these three parameters.

IV-E Regression Performance Comparison

The average RMSEs and CCs for the five algorithms across the 15 subjects are shown in Fig. 4, and the RMSEs and CCs for the individual subjects are shown in Fig. 5. Observe that DAMF, OwARR and OwARR-SDS had very similar CCs, and all of them were better than the CCs of BL1 and BL2. Additionally:

Except for BL1, whose model does not depend on $m$, all the other four algorithms gave smaller RMSEs as $m$ increased, which is intuitive. 2. 2.

BL1 had the smallest RMSE when $m=0$. However, as $m$ increased, DAMF, OwARR and OwARR-SDS quickly outperformed BL1. This suggests that there is large individual difference among the subjects, and hence a subject-independent model is not desirable. 3. 3.

Because BL2 used only subject-specific calibration data, it cannot build a model when $m=0$, i.e., when there was no subject-specific calibration data at all. However, all the other four methods can, because they can make use of data from other subjects. BL2’s performance was the worst, because it cannot get enough training when there is only a small number of subject-specific calibration samples. 4. 4.

OwARR and OwARR-SDS had almost identical average RMSEs, which were smaller than those of BL1, BL2 and DAMF. More importantly, the RMSEs of OwARR and OwARR-SDS almost converged as soon as the first batch of subject-specific samples were added, suggesting that they only need very few subject-specific samples to train, which is very desirable in practical calibration.

In summary, the three DA based approaches generally had better performance than BL1, which does not use subject-specific data at all, and also BL2, which does not use auxiliary data at all. This suggests that DA is indeed beneficial. Moreover, our proposed OwARR and OwARR-SDS achieved the best overall performances among the five algorithms.

We also performed two-way Analysis of Variance (ANOVA) for each $m$ to check if the RMSE differences among the five algorithms were statistically significant, by setting the subjects as a random effect. Two-way ANOVA showed statistically significant differences among them ($p<0.01$) for all $m$. Then, non-parametric multiple comparison tests using Dunn’s procedure [8, 9] 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 [2]. The $p$-values are shown in Table I, where the statistically significant ones are marked in bold. Observe that the differences between OwARR and the other three algorithms (BL1, BL2, and DAMF) were always statistically significant when $m>0$, so were the differences between OwARR-SDS and the two baseline algorithms. The difference between OwARR-SDS and DAMF was statistically significant for $m\in[5,75]$. There was no statistically significant difference between OwARR and OwARR-SDS.

Finally, we can conclude that given the same amount of subject-specific calibration data, OwARR and OwARR-SDS can achieve significantly better estimation performance than the other three approaches. Or, in other words, given a desired RMSE, OwARR and OwARR-SDS require significantly less subject-specific calibration data than the other three approaches. For example, in Fig. 4, the average RMSE for BL2 when $m=100$ was $0.2988$, whereas OwARR and OwARR-SDS can achieve even smaller RMSEs without using any subject-specific calibration samples. The average RMSEs for OwARR and OwARR-SDS when $m=5$ were $0.2347$ and $0.2348$, respectively, whereas DAMF needed at least 45 subject-specific samples to achieve these RMSEs, and BL2 needed at least 100 samples.

IV-F Computational Cost

In this subsection we compare the computational cost of the five algorithms, particularly, OwARR and OwARR-SDS, because the primary goal of SDS is to down-select the number of auxiliary subjects and hence to reduce the computational cost of OwARR.

Fig. 6 shows the average number of similar subjects selected by SDS for the 15 subjects. Observe that most of the time fewer than seven subjects (half of the number of auxiliary subjects) were selected.

To quantify the computational cost of the five algorithms, we show in Fig. 7 the training times for different $m$, averaged over 10 runs and across the 15 subjects. The platform was a Dell XPS 13 notebook, with Intel Core i7-5500M CPU@2.40GHz, 8GB memory, and 256GB solid state drive. The software was Matlab R2015b running in 64-bit Windows 10 Pro. Each algorithm was optimized to the best ability of the authors. Observe that the training time of BL1, BL2, and DAMF was almost constant, whereas the training time of OwARR increased monotonically as $m$ increased. Interestingly, the training time of OwARR-SDS decreased slightly as $m$ increased, because Fig. 6 shows that generally the average number of similar subjects selected by SDS decreased as $m$ increased.

The computational costs of OwARR and OwARR-SDS were much higher than DAMF, because they used more sophisticated DA approaches. However, except for $m=0$, at which point OwARR and OwARR-SDS had identical training time, the training time of OwARR-SDS was on average only about 49% of OwARR. This 51% computation time saving is very worthwhile when the number of source domains is very large and hence computing OwARR for all the source domains is too slow.

We also investigated the scalability of OwARR with respect to $Z$, the number of source domains, and $n$, the number of samples in each source domain. Because we only had 14 source domains in this dataset, we bootstrapped them to create additional domains when $Z\geq 14$. The results are shown in Fig. 8. Observe from Fig. 8 that the computational cost of OwARR increased linearly with the number of source domains, which is intuitive, because OwARR performs DA for each source domain separately and then aggregates the results. However, Fig. 8 shows that the computational cost of OwARR increased superlinearly with the number of samples in the source domains. Least-squares curve fitting found that the computation time was about $0.0000021\cdot n^{1.8}+0.035$ seconds, i.e., the computational cost is $O(n^{1.8})$ for 14 source domains.

Finally, it is important to note that the above analyses are only for the training of the algorithms. Once the training is done, the resulting OwARR and OwARR-SDS models can be executed much faster.

IV-G Robustness to Noises

It is also important to study the robustness of the five algorithms to noises. According to [68], 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 [68], for each model input, we randomly replaced $q\%$ ($q=0,\,10,...,50$) of all epochs from the new subject with a uniform noise between its minimum and maximum values. After this was done for both the training and testing data, we trained the five algorithms on the corrupted training data and then tested their performances on the corrupted testing data. The RMSEs for three different $m$ (the number of labeled subject-specific samples), averaged across 15 subjects with five runs per subject, are shown in Fig. 9. Observe that as the noise level $q$ increased, generally all algorithms had worse RMSEs. OwARR and OwARR-SDS still had the smallest RMSEs among the five when $q$ was small. However, when $q$ increased, DAMF became the best. This suggests that OwARR and OwARR-SDS may not be as robust as DAMF with respect to attribute noises, but when the noise level is low, the performance improvement achieved from the sophisticated optimizations in OwARR and OwARR-SDS dominates, and hence they are still the best algorithms among the five. When the noise level is high, we may need some noise handling approaches, e.g., noise correction [68], before applying OwARR and OwARR-SDS.

IV-H Parameter Sensitivity Analysis

The OwARR algorithm has three adjustable parameters: $\sigma$, which determines the weight $w^{t}$ for the target domain samples; $\lambda$, which is a regularization parameter minimizing the distances between the marginal and conditional probability distributions in the source and target domains; and $\gamma$, which maximizes the approximate Pearson correlation coefficient between the true and estimated outputs. It is interesting to study whether all of them are necessary.

For this purpose, we constructed three modified versions of the OwARR algorithms by setting $\sigma$, $\lambda$ and $\gamma$ to zero, respectively, and compared their average RMSEs with that of the original OwARR. The results are shown in Fig. 10. Observe that the original OwARR had better performance than all three modified versions, suggesting that all three parameters in OwARR contributed to its superior performance.

Next we studied the sensitivity of OwARR to the three adjustable parameters, $\sigma$, $\lambda$ and $\gamma$. The results are shown in Fig. 11(a)-(c). Observe that OwARR is robust to $\sigma$ in the range of $[0.1,0.4]$, to $\lambda$ in the range of $[1,20]$, and to $\gamma$ in the range of $[0.01,1]$.

Additionally, three type-1 triangular FSs have been used in conditional probability distribution adaptation (Section II-E) in this paper for simplicity. It is also interesting to study the sensitivity of the OwARR algorithm to the number of FSs. The results are shown in Fig. 11. Observe that OwARR gives the optimal performance when the number of FSs is between 2 and 5, but its performance gradually deteriorates when the number of FSs further increases. This is intuitive, because the target domain has a limited number of labeled training samples, so as the number of FSs increases, the number of target domain samples that fall into each fuzzy class decreases, and hence the computed fuzzy class means are less reliable. As a result, the distance between the conditional probability distributions [see (20)] cannot be reliably computed.

Another interesting questions is: what would be the performance of OwARR if no FSs are used at all, i.e., conditional probability distribution adaptation is disabled? This is corresponding to the left-most slice in Fig. 11, where the number of FSs is zero. Observe that this results in worse RMSEs than the case that two to five FSs are used in conditional probability distribution adaptation, suggesting the FS approach is beneficial.

IV-I Effectiveness of the Ensemble Fusion Strategy

From Algorithm 1 (Algorithm 2) it is clear that the final step of OwARR (OwARR-SDS) uses ensemble learning: the base DA models are aggregated using a weighted average to obtain the final regression model, and the weight is inversely proportional to the training RMSE of the corresponding base DA model. In this subsection we study whether this fusion strategy is effective. The performances of the 14 base DA models and the final aggregation model for a typical subject are shown in Fig. 12. Observe that the aggregated model is better than most base DA models, and is also close to the best base DA model (which is unknown in practice), suggesting that the fusion strategy is effective. However, it may be possible that a better fusion strategy can make the final model outperform all base DA models. This will be one of our future research directions.

V Conclusions and Future Research

Transfer learning, which improves learning performance in a new task by leveraging data or knowledge from other relevant tasks, represents a promising solution for handling individual differences in BCI. Previously we have proposed a weighted adaptation regularization (wAR) algorithm [61, 59] for offline BCI classification problems, an online weighted adaptation regularization (OwAR) algorithm [60] for online BCI classification problems, and a SDS approach [60, 61] to reduce the computational cost of wAR and OwAR. In this paper we have proposed an OwARR algorithm to extend the OwAR algorithm from classification to regression, and validated its performance on online estimation of driver drowsiness from EEG signals. Meanwhile, we have also extended the SDS algorithm for classification in [60] to regression problems, and verified that OwARR-SDS can achieve similar performance to OwARR, but save about half of the computation time. Both OwARR and OwARR-SDS use fuzzy sets to perform part of the adaptation regularization, and OwARR-SDS also uses fuzzy sets to select the closest source domains.

Though OwARR and OwARR-SDS have demonstrated outstanding performance, they can be enhanced in a number of ways, which will be considered in our future research. First, Fig. 5 shows that OwARR and OwARR-SDS had worse RMSEs than BL1 for some subjects. This indicates that they still have room for improvement: we could develop a mechanism to switch between BL1 and OwARR (OwARR-SDS) so that a more appropriate method is chosen according to the characteristics of the new subject, similar to the idea of selective TL [51]. Second, we will extend OwARR and OwARR-SDS to offline calibration, where the goal is to automatically label some initially unlabeled subject-specific samples with a small number of queries [61]. Semi-supervised learning can be used here to enhance the learning performance. Third, in this paper we combine the base learners using a simple weighted average, where the weights of the base learners are inversely proportional to their corresponding training RMSEs. This may not be optimal because what really matter here are the testing RMSEs. In online calibration it is not easy to estimate the testing RMSEs because we do not know what samples will be encountered in the future; however, in offline calibration we can better estimate the testing performances of the base learners using a spectral meta-learner approach [58], and hence a better model fusion strategy could be developed. Fourth, similar to offline classification problems [55, 59, 29], in offline regression problems we can also integrate DA with active learning [40, 57] to further reduce the offline calibration effort. Finally, we will apply the online and offline DA algorithms to other regression problems in BCI and beyond to cope with individual differences, e.g., estimating the continuous values of arousal, valence and dominance from speech signals [64] in affective computing.

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.

[Fuzzy Sets (FSs)]

FS theory was first introduced by Zadeh [67] in 1965 and has been successfully used in many areas, including modeling and control [49, 65], data mining [35, 66, 62], time-series prediction [46, 14], decision making [30, 31, 37], etc.

A FS $X$ is comprised of a universe of discourse $D_{X}$ of real numbers together with a membership function (MF) $\mu_{{}_{X}}:\ D_{X}\to[0,1]$, i.e.,

[TABLE]

Here $\int$ denotes the collection of all points $x\in D_{X}$ with associated membership degree $\mu_{{}_{X}}(x)$. An example of a FS is shown in Fig. 13. The membership degrees are $\mu_{X}(1)=0$, $\mu_{X}(3)=0.5$, $\mu_{X}(5)=1$, $\mu_{X}(6)=0.8$, and $\mu_{X}(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).

FSs are frequently used in modeling concepts in natural language, which may not have clear boundary. For example, we may define a hot day as temperature equal to or above $30^{\circ}$C, but is $29^{\circ}$C hot? If we represent hot as a binary set $\{x|x\geq 30\}$, then $29^{\circ}$C is not hot, because it does not belong to the binary set hot. However, this does not completely agree with people’s intuition: $29^{\circ}$C is very close to $30^{\circ}$C, and hence it is somewhat hot. If we represent hot as a FS, we may say $29^{\circ}$C is hot with a membership degree of $0.9$, which sounds more reasonable.