MENTAL FATIGUE MONITORING USING BRAIN DYNAMICS PREFERENCES
Year: 2018

Abstract:
Driver’s cognitive state of mental fatigue significantly affects driving performance and more importantly public safety. Previous studies leverage the response time (RT) as the metric for mental fatigue and aim at estimating the exact value of RT using electroencephalogram (EEG) signals within a regression model. However, due to the easily corrupted EEG signals and also non-smooth RTs during data collection, regular regression methods generally suffer from poor generalization performance. Considering that human response time is the reflection of brain dynamics preference rather than a single value, a novel model called Brain Dynamics ranking (BDrank) has been proposed. BDrank could learn from brain dynamics preferences using EEG data robustly and preserve the ordering corresponding to RTs. BDrank model is based on the regularized alternative ordinal classification comparing to regular regression based practices. Furthermore, a transition matrix is introduced to characterize the reliability of each channel used in EEG data, which helps in learning brain dynamics preferences only from informative EEG channels. In order to handle large-scale EEG signals and obtain higher generalization, an online-generalized Expectation Maximum (OnlineGEM) algorithm also has been proposed to update BDrank in an online fashion. Comprehensive empirical analysis on EEG signals from 44 participants shows that BDrank together with OnlineGEM achieves substantial improvements in reliability while simultaneously detecting possible less informative and noisy EEG channels.

heading: 1 INTRODUCTION
text: As reported by sleep health report (Adams et al., 2017), mental fatigue is a major cause in 33%− 45% of all road accidents. In general, mental fatigue (Boksem & Tops, 2008) refers to the inability to maintain optimal cognitive performance in continuous task of the high demand of cognitive activity. Such inability in the context of driver could lead to accidents with severe consequences. Individuals may find themselves in a mental fatigue state because of lack of sleep, continuous driving for longtime, midnight driving, monotonous driving, and driving during the influence of sleeping drugs or sleep disorders (Ji et al., 2004; Ting et al., 2008).
In response to these critical issues, several methods (Jap et al., 2009; Wascher et al., 2014; Cook et al., 2007; Lal et al., 2003; de Naurois et al., 2017) have been proposed to estimate and predict the mental fatigue based on EEG and RT (Fig. 1a). Some of these methods, however, performed considerably well for some participants but failed for others due to lack of robustness. There are several challenges behind such instability and one of such problem is how to use RT effectively. RT is the most resourceful piece of information to predict mental fatigue. However, it is easily affected by the instrumental error, mind wandering or any other task non-related factor. A previous study (Wei et al., 2015) tried to overcome this problem by providing different techniques to adjust RT, they also tried removing outliers nevertheless, failed to make it work for all participants. The regression assumption of this method between EEG signals and RT is not correct. Human’s RT is usually the result of preference (Izuma & Adolphs, 2013) in brain dynamics during the task, than just a single value. These preferences of human can be affected by different cognition (Möckel et al., 2015) like mind wandering (Lin et al., 2016), and/or lower level of attention (Chuang et al., 2018). Therefore, the relationship between EEG and RT including the extreme/abnormal RT should be taken care in the way that reflect human brain dynamics preferences by the developed technique itself.
Another important problem is related to different brain regions, which are normally responsible for different functionalities. There was an attempt to choose different brain regions (Wascher et al., 2014) for a method during mental fatigue prediction but these regions of the brain are not necessary same for all participants (Gramann et al., 2006). For example, existing work (Wascher et al., 2014) used frontal theta to represent a different level of mental fatigue for all participants. In such case learning model’s reliability would inevitably degrade because of possibly noisy channels chosen, on different brain regions, by the method. Some previous work (de Naurois et al., 2017), attempted to solve this issue using artificial neural network models but still failed to provide convincing results. Again these mentioned work reflect developed methods should be based on brain dynamics preferences rather than fix model or regions of brain.
To overcome the above-mentioned problems, a new approach has been proposed. We called it BDrank such that brain dynamics ranking. This approach not only learns from brain dynamics preferences for mental fatigue but also other cognitive states (Lal et al., 2003), while effectively preserving the exact ordering of RT (Fig. 1b). This approach surprisingly improved over defects of previous models and their performance due to noisy and extreme RT. Furthermore, BDrank model also proposes to use transition matrix to evaluate the high confidence BDrank (HC-BDrank) sources among different EEG channels, which contributes highly toward task performance. In order to handle large-scale EEG signals and obtain higher generalization, an online-generalized Expectation Maximum (OnlineGEM) (Roche, 2011) algorithm also has been proposed. Comprehensive empirical experiments on EEG signals from 44 participants show that BDrank model together with OnlineGEM algorithm delivers substantial improvements and robust performance while simultaneously providing the information about noisy or less informative EEG channels.

heading: 2 MATERIALS AND METHODS:
text: Huang et al. Information transfer during drowsy driving
FIGURE 1 | Sustained-attention driving task. (A) Event-related lane-deviation paradigm. (B–C) EEG and behavior were recorded simultaneously. (D) The
experimental paradigm was implemented in a VR-based driving simulator, where a 3-D surrounding view simulated a monotonous highway scene and a vehicle was
mounted on a six degree-of-freedom dynamic platform. (E) The car speed was fixed at 100 km/h and the subject was asked to keep the car cruising on the central of
the lane.
processes X and Y. We assume that each of these processes can be approximated by a stationary Markov Process of finite order d. Thus, we can reconstruct the state space of the process X by a delay embedded vector of dimension d with past values. Representing the two time series as X = xt and Y = yt, the delay embedded vector is defined as x d t = (
xt, xt−τ , xt−2τ , · · · , xt−(d−1)τ ) ; similar representation can also
be made for ydt . The dimension of the embedding space is d, and the delay is τ . Under the assumption that the system X can be approximated by a stationary Markov process of order d, the transition probabilities that describe the system are given by:
p(xt+1
∣ ∣ ∣ xdt ). (1)
The entropy rate of the systemX is the average number of bits that is required to represent an additional state provided all previous states are known. Thus, the entropy rate can be computed as follows
h(xt+ u
∣ ∣ ∣ xdt ) = − ∑
xt+ u,x d t
p(xt+ u, x d t ) log p(xt+u
∣ ∣ ∣ xdt ), (2)
where p(xt+u
∣ ∣ ∣ xdt ) = p(xt+u, x d t )/p(x d t ),
u is the prediction time, and p(∗) is the probability. If the two processes are independent, there will be no transfer
of information and p(xt+u
∣ ∣ ∣ xdt ) = p(xt+u ∣ ∣ ∣ xdt , y m t ). As proposed
by Schreiber (2000), a measure of deviation from this generalized Markov property can be computed using Kullback divergence or mutual information and that is a directedmeasure of information transfer from Y to X. The amount of information transferred (i.e., transfer entropy) from process Y to process X, denoted
as TE(Y → X), is computed as (Schreiber, 2000; HlaváèkováSchindler et al., 2007; Vicente et al., 2011)
TE(Y → X) = ∑
p(xt+u, x d t , y m t ) log
p(xt+u
∣ ∣ ∣ xdt , y m t )
p(xt+u
∣ ∣ ∣ xdt )
, (3)
where ymt = { yt, yt−τ, · · · , yt−(m−1)τ }
indicating that the process X depends onm states of Y. The Equation (3) can be rewritten in terms of differential entropy as the following
TE(Y → X) = H(xdt , y m t )−H(xt+u, x d t , y m t )
+ H(xt+u, x d t )−H(X d t ). (4)
Note that TE is inherently asymmetric, that is TE(Y → X) 6= TE(X → Y). In addition, when the processes are mutually independent then TE(Y → X) = TE(X → Y) = 0.
The estimation of parameters and the calculation of TE were performed by TRENTOOL (version 2.0.4) (Lindner et al., 2011). Specifically, the TE values were estimated by the k-nearest neighbor approach (Kraskov et al., 2004; Vicente et al., 2011), with k = 4 as suggested by Kraskov et al. (2004). The embedding delay (τ) was determined based on the Cao criterion (Cao, 1997), and the dimension (d andm) was obtained by an effective search algorithm (Cormen et al., 2001; Lindner et al., 2011; Vicente et al., 2011). Theiler correction tries to remove autocorrelation effect from the density estimation. For the nearest neighbor search, it discards all samples which are very close in time with respect to a reference point. Here the Theiler correction window (T) was set to 1. In this study the prediction time u was set to 5. This was determined by finding the maximum value of TE over a set of choices: (5, 10, 20, 40, 60, 80, and 100ms).
Frontiers in Human Neuroscience | www.frontiersin.org 3 October 2015 | Volume 9 | Article 570
Figure 2: Sustained-attention riving tas
a) Experiment Paradigm: This paper utilized the 33-channels EEG data recorded in the previous study (Huang et al., 2015) from 44 adult participants while performing long sustained attention task. The experiment has been conducted using a virtual-reality (VR) dynamic driving simulator (Fig. 2D-E). The task involves driving on the fourlane highway while lane-departure events were randomly induced deviation toward the side of the road from the original position. Each participant was instructed to quickly respond to steer back to original position. A complete trial
in this study (Fig. 2A), includes 10s baseline, deviation onset, response onset, and response offset (Fig. 2B-C). The next trial occurs within an interval of 5-10s after finishing the current trial. Each participant completed T trials within 1.5h. For each trial i, the EEG signals {xn,i}Nn=1 from N different channels were recorded simultaneously and the corresponding response time RTi was also collected afterwards. If a participant fell asleep during the experiment, there was no feedback to wake him up.
2
The response time is an intuitive indicator used to assess human mental fatigue. Therefore, the common practice for mental fatigue monitoring is to find a robust mapping for humans’ response time (RT) to an emergent situation using the EEG signals recorded beforehand. The natural way to forecast the response time with the EEG signals is to formulate it as a regression task (Fig. 1a), namely finding a nonlinear mapping (e.g. neural network, SVR) from the EEG signals x to the corresponding RT. However, due to the easily corrupted properties of the EEG signals and the existence of extreme values in response times during data collection, focusing on predicting the exact value of a noisy and non-smooth measurement (the response time) is easier to create a near-perfectly fitted model with poor generalization performance (See Table 1 and Fig. 3 for more details). This creates a dilemma: it requires a powerful learning model to predict response time with the complex EEG signals (indeed, it is exactly our target) but it is not so significant to excessively approximate the exact value of response time, especially the extreme values.
Here comes the problem, how to find an efficient way to learn from the noisy response time while the exact value is not so significant? Actually, the RTs are defined in the totally ordered space R. The totally ordered space owns its structure meanings, which are preserved by the pairwise comparisons between the RTs. The pairwise comparisons indeed preserve the whole relative structure information between the response times while ignoring their absolute numerical information. Therefore, predicting the orderings of the pairwise comparisons can be accepted as a regularized alternative of previous regression model. Further, it makes it more flexible to consider more powerful learning unit in the proposed model.
Therefore, there is no longer requirement to estimate the RT with a regression model. Instead, the proposed approach will transform it into an ordinal classification problem, and focus on correctly preserving the whole orderings between the pairwise RT comparisons (Fig. 1b). First, the preference propositions1 could be constructed as follows,
D1 = {ρm}M1m=1 = {RTm,1 > RTm,2} M1 m=1, D2 = {ρm′} M2 m′=1 = {RTm′,1 ≈ RTm′,2} M2 m′=1. (1)
where M1 and M2 denote the number of type-1 and type-2 preference propositions. D1 denotes the type-1 preference propositions that the orderings between the RTs are significant; D2 denotes the type-2 preference propositions that the RTs in each comparison are comparable. Then, the brain dynamics preference was constructed for each proposition with the corresponding pairwise EEG signals recorded from each channel. For brevity of notations, new index (the notation in Eq.1) is adopted in the following, instead of the original index used in sustained-attention driving task (Fig. 1). Namely, the m-th proposition ρm ∈D1 ⋃ D2 denotes the pairwise comparison RTm,1 > RTm,2 or RTm,1 ≈ RTm,2. And the pairwise brain dynamics preference (x1n,m,x2n,m) denotes the features recorded within the n-th channel for each preference proposition ρm ∀m = 1,2, · · · ,M1 +M2. In this paper, the 10s baseline (Fig. 2B) as the feature vector has been adopted, which is assumed to be long enough to detect any significant changes in brain activity (Zhang, 2000). This followed by exploring the relation between the 10s baseline x (∈ Rk) and the preference proposition ρm under the following four assumptions: (1) different participants are independent during the data collection process; (2) Different EEG sensors used for recording are recorded independently from scalp without influencing other sensors (Homan et al., 1987; Teplan et al., 2002); (3) Different trials are conducted independently during the data collection process; (4) The collected response time are slightly corrupted by inherent (basically irremovable) sources of noise, but the ranking relationships are preserved to some extents.
b) Brain Dynamics Preferences: Again, instead of modelling the dependence between the RT and the 10s baseline as a regression problem, we aim to predict the orderings of the pairwise comparisons using the brain dynamics preferences, namely f (x1,x2)→ ρ . Further, a preference proposition ρ has three states: 1,0,−1, denoting win (RT1 > RT2), tie (RT1 ≈ RT2) and loss (RT1 < RT2), respectively. However, classical classification models, e.g. logistic ordinal regression (Harrell, 2001), are infeasible for our problem, due to the lack of a normalized probability definition for three states. Since two types of preference propositions D1 and D2 are considered in our problem, we tailor-define a normalized probability definition, namely first normalizing the probabilities of the states (1,−1)
1We used the term “preference” intentionally to show that brain dynamics keep changing w.r.t. human behaviours and it happens because the human brain prefers one decision over others. Therefore, we prefer to call it“preference” than “classification”.
(exclusively to D1) to 1, then generalize the probability definition to state 0. To be specific, it can be mathematically formulated as follows,
P(ρ|w,x1,x2) =  σ(wT ∆x)[1−κ(wT ∆x)] ρ = 1 κ(wT ∆x) ρ = 0 σ(−wT ∆x)[1−κ(wT ∆x)] ρ =−1
(2)
where σ(z) = 1/(1+ e−z) is the sigmoid function and σ(−z) = 1−σ(z). The symbol ∆x denotes the subtraction (x1− x2) between the brain dynamics preferences (x1,x2). Following Weng & Lin (2011), the probability of a tie is modelled as the geometric mean between a win and a loss, namely κ(wT ∆x) = √ σ(wT ∆x)[1−σ(wT ∆x)].
c) Preference State Transition: Furthermore, considering the different functions of different regions in the human brain, the relative contributions of different channels to human response time may vary a lot. For example, the information conveyed by positive channels is positive related to the RT, while negative channels may convey the information which is negative related to the RT. There are also some noisy (non-relevant) channels which are independent to the learning task. Therefore, if we directly model the EEG preferences recorded in each channel without making any distinctions about the channel state (i.e. positive, noisy and negative), the model’s reliability would inevitably degrade. Note that a channel is called as “noise” if the current algorithms could not extract useful brain information with EEG signals from this channel (Alharbi, 2018; Lin et al., 2018).
In the following, a transition matrix Πn is introduced to characterize the reliability of each channel n w.r.t the corresponding task. Let ρ denote the preference proposition and ρ(n) denote the prediction from the n-th channel. They are defined on the finite state space S = {1,0,−1}. Then we have
Πn = P(ρ|ρ(n)) =  πn 0 (1−πn)0 1 0 (1−πn) 0 πn  , (3) where Pi, j(ρ|ρ(n)) =P(ρ = S j|ρ(n) = Si). Note that we do not consider the transition between the type1 and type-2 preference propositions (i.e. P(ρ = {1,−1}|ρ(n) = 0) = 0 and P(ρ = 0|ρ(n) = {1,−1}) = 0), since the equal cases are hard to measure when do prediction. A promising approach to generalize the transition matrix Πn (Eq.3) is to introduce the concept of the confidence region to measure the equal cases (Pregibon et al., 1981).
Remark: The parameter πn in the transition matrix Πn actually indicates the reliability of the n-th channel ∀n = 1,2, · · · ,N. It additionally helps to divide the channels into three states: (1) positive channels with πn close to 1, the ranking model (Eq.2) can extract enough information from the n-th channel, and exactly predict the state of the preference proposition. (2) Noisy channels with πn approximating to 0.5, the ranking model can not extract any useful information from the n-th channel. (3) Negative channels with πn close to 0, the ranking model can extract enough information from the n-th channel, but the prediction states are exactly opposite to the proposition states. The identified positive and negative channels are all considered as informative EEG channels, which helps in learning reliable models for the corresponding task.
With the introduced transition matrix Πn, the marginal likelihood function for each proposition ρ can be further represented as P(ρ) = Eρ(n) [ P(ρ|ρ(n))P(ρ(n)) ] . Specifically,
P(ρ|w,Πn,x1,x2) =  [πnσ(wT ∆x)+(1−πn)σ(−wT ∆x)][1−κ(wT ∆xn,m)] ρ = 1
κ(wT ∆xn,m) ρ = 0 [(1−πn)σ(wT ∆x)+πnσ(−wT ∆x)][1−κ(wT ∆xn,m)] ρ =−1
(4)
For sake of simplicity, the subscripts (n,m) are omitted. Let D =D1 ⋃
D2 denotes all the preference propositions and X represents the recorded EEG signals from N different channels. We further extend our robust ordinal classification (Eq.4) for BDrank into the Bayesian version. A Gaussian prior is introduced for w (i.e., w ∼ N(µ,Σ)). Since the transition matrix Πn only depends on the parameter πn, we only focus on estimate the parameter πn ∀n = 1,2, · · · ,N in the following. Let π denote {πn}Nn=1, and we introduce a Beta prior for each πn (i.e., π ∼ B(α,β ) = ∏Nn=1 B(αn,βn)). Then,
our BDrank can be formulated into a maximum a posteriori (MAP) estimate as follows,
argmax π,w P0(π)P0(w)P(D |π,w,X) = P0(π)P0(w)P(D1|π,w,X)P(D2|π,w,X)
= P0(π)P0(w) N
∏ n=1 [ M1 ∏ m=1 P(ρm = 1|πn,w,∆xn,m) M2 ∏ m′=1 P(ρm′ = 0|πn,w,∆xn,m′) ]
(5)
= B(π|α,β )N(w|µ,Σ) N
∏ n=1 [ M1 ∏ m=1 ( [πnσ(wT ∆x)+(1−πn)σ(−wT ∆x)][1−κ(wT ∆xn,m)] ) M2 ∏ m′=1 κ(wT ∆xn,m′) ] .
Note that due to the symmetry of the state probability (Eq.2) and transition matrix (Eq.3) w.r.t. state 1 and −1, the resultant marginal likelihood (Eq.4) remains symmetry w.r.t. state 1 and −1. For simplicity, D1 is constructed using the preference propositions with state ρ = 1 only. M1 and M2 denote the number of type-1 and type-2 preference propositions, namely |D1|= M1 and |D2|= M2. The variable n iterates over the channels. m and m′ iterate over two types of preference propositions, respectively. Now our aim is to estimate the model parameters (w and π) by maximizing Eq.5. In principle, any solution strategies for MAP can be considered to solve this problem. See Section 3 for optimization details.
d) Reliability analysis and channel state estimation Note that our BDrank (Eq.5) indeed trains a mixture of two complementary classifiers, which share the same parameter w. It is different from classical mixture models, since it clusters at the channel level instead of the sample level. In particular, in terms of the positive channels with πn close to 1, BDrank relies on the first classifier to update the shared parameter w. In terms of the negative channels with πn close to 0, Eq.5 automatically switches to the opposite classifier which can extract correct information from the negative channels and update the shared parameter w accordingly. Further, BDrank is robust to the noisy channels with πn approximating to 0.5, because Eq.5 gives up extracting information from the noisy channels by assigning a constant likelihood (i.e. 0.5) to each brain dynamics preference. Note that the estimated πn can be leveraged as the indicator to detect noisy channels with πn ≈ 0.5,∀n = 1,2, · · · ,N. See Fig. 4 in the experimental section for more details.

heading: 3 OPTIMIZATION METHODS
text: In this section, we describe a generalized Expectation-Maximization (EM) algorithm (Dempster et al., 1977) to solve the proposed BDrank (Eq.5). Since the feasible region of πn is restricted to [0,1], the gradient-based optimization methods would make our solution inaccurate and inefficient. The EM algorithm is an efficient iterative procedure to compute the MAP problem in presence of latent variable (ρ(n)m in Eq.5). EM avoids calculating the derivative to the expectation of the latent variable directly, and resorts to a surrogate lower bound to optimize. Therefore, EM, a silver bullet for MAP with latent variable, can significantly simplify the optimization over parameter πn for Eq.5.

heading: 3.1 GENERALIZED EM FOR BDRANK
text: For each type-1 proposition ρm, we introduce an auxiliary variable δ (n) m (∈ {1,0}) for each channel, representing the consistency between the preference proposition ρm and the prediction ρ (n) m given by the n-th channel. Specifically, δ (n)m = 1 denotes the prediction ρ (n) m given by the first classifier is consistent with the preference proposition ρm, and δ (n) m = 0 denotes the prediction ρ (n) m estimated by the second classifier is consistent with the preference proposition ρm. We can therefore find an equivalent formulation of our BDrank model involving the auxiliary variable Ξ = {δ (n)m }M1m=1.
P(D ,Ξ,π,w|X) = P0(π)P0(w)P(D1,Ξ|π,w,X)P(D2|w,π,X)
= P0(π)P0(w) N
∏ n=1 [ M1 ∏ m=1 P(ρm = 1,δ (n) m |πn,w,∆xn,m) M2 ∏ m′=1 P(ρm′ = 0|πn,w,∆xn,m′) ]
(6)
= P0(π)P0(w) N
∏ n=1 [ M1 ∏ m=1 [ πnσ(wT ∆xn,m) ]δ (n)m [(1−πn)σ(−wT ∆xn,m)]1−δ (n)m [1−κ(wT ∆xn,m)] M2∏ m′=1 κ(wT ∆xn,m′) ] .
Now, we can deal with the joint distribution directly, which leads to significant simplifications for optimization. The complete log likelihood can be written as
logP(D ,Ξ,π,w|X) = logP0(π)+ logP0(w)+ N
∑ n=1
M1
∑ m=1
log[1−κ(wT ∆xn,m)]+ N
∑ n=1
M2
∑ m′=1 logκ(wT ∆xn,m′) (7)
+ N
∑ n=1
M1
∑ m=1
[ δ (n)m logπnσ(wT ∆xn,m)+(1−δ (n) m ) log(1−πn)σ(−wT ∆xn,m) ] .
Expectation Step In the expectation step, we first calculate the expected value of the auxiliary variable δ (n)m w.r.t its posterior distribution P(δ (n)m |π,w,ρm = 1,xn,m) ∀n = 1,2, · · · ,N,∀m = 1,2, · · · ,M1:
E P(δ (n)m |π,w,ρm=1,xn,m)
[δ (n)m ] = P(ρm = 1,δ
(n) m |πn,w,∆xn,m)
P(ρm = 1|πn,w,∆xn,m) =
1
1+ (1−πn)σ(−w T ∆xn,m)[1−κ(wT ∆xn,m)]
πnσ(wT ∆xn,m)[1−κ(wT ∆xn,m)]
, γ(δ (n)m ),
(8) where γ(δ (n)m ) denotes the degree of the consistency between the prediction ρ (n) m and the preference
proposition ρm. Then, the expected value of the complete-data log likelihood function w.r.t the posterior expectation of the auxiliary variable Ξ = {δ (n)m }M1m=1 can be represented as follows:
L (w,π) = N
∑ n=1
[(αn−1) logπn +(βn−1) log(1−πn)]− 1 2 (w−µ)T Σ−1(w−µ)+
N
∑ n=1
M1
∑ m=1 log[1−κ(wT ∆xn,m)]
+ N
∑ n=1
M2
∑ m′=1
logκ(wT ∆xn,m′)+ N
∑ n=1
M1
∑ m=1
[ γ(δ (n)m ) logπn +(1− γ(δ (n) m )) log(1−πn) ] (9)
+ N
∑ n=1
M1
∑ m=1
[ γ(δ (n)m ) logσ(wT ∆xn,m)+(1− γ(δ (n) m )) logσ(−wT ∆xn,m) ] + constant.
Maximization Step In the maximization step, we maximize the objective function Eq.9 w.r.t the model parameters π and w, respectively. In terms of π , we set the gradient of Eq.9 w.r.t πn to zero and obtain the following estimate for πn:
πnewn = ∑M1m=1 γ(δ (n) m )+αn−1
M1 +αn +βn−2 , n = 1,2, · · · ,N. (10)
In terms of w, due to the complexity of the sigmoid function, we do not have a closed form solution for w and we need to use gradient-based optimization methods. In the following, we adopt the LBFGS to optimize w, the objective function L (w) consists of the second, the third, the fourth and the sixth parts of Eq.9 and the gradient function g(w) can be represented as follows,
g(w) =−Σ−1(w−µ)+ N
∑ n=1
M1
∑ m=1 1−2σ(wT ∆xn,m) 2(1− 1κ(wT ∆xn,m) ) ∆xn,m + N ∑ n=1 M1 ∑ m=1 (γ(δ (n)m )−σ(wT ∆xn,m))∆xn,m
+ 1 2
N
∑ n=1
M2
∑ m′=1 (1−2σ(wT ∆xn,m′))∆xn,m′ . (11)
wnew can be solved with L-BFGS using L (w) and g(w). The EM algorithm then iterates the E-step and M-step until convergence is achieved.

heading: 3.2 ONLINE GEM FOR BDRANK
text: The generalized EM approach introduced in Section 3.1 is inefficient for large-scale datasets, because we need to iteratively calculate the gradient with respect to parameters π and w over all samples during each maximization step. Motivated from the stochastic approximation literature, we introduce an online-generalized Expectation Maximization (OnlineGEM) approach, which resorts to stochastic mini-batch optimisation to learn the parameters. To be specific, OnlineGEM approximates the updated π and w in batch EM with a single sample or mini-batch samples. Since a mini-batch samples cannot be a perfect approximation for the whole dataset, we interpolate between the new and former estimations with a decreasing step-size ηk2, as in Liang & Klein (2009).
2ηk = (k+ 2)−τ0 , where k is the number of iterations and 0.5 < τ0 < 1. The smaller the τ0, the larger the update ηk, and the more quickly we forget (decay) our old parameters. This can lead to swift progress but also generates instability.
Before the k-th iteration, we randomly downsample a mini-batch Dk from the preference propositions D . The number of two types preference propositions in Dk are M∗1 and M ∗ 2 , which are much smaller than the corresponding total size M1 and M2, respectively.
The expectation step remains similar. The difference is that we only need to calculate the posterior expectation of the auxiliary variable δ (n)m over the mini-batch Dk.
In the maximization step, we maximize the objective function, calculated on the mini-batch Dk, with regard to model parameters π and w. In terms of parameter πn, since its marginal distribution belongs to exponential family, we could perform the stochastic update in the space of sufficient statistics (Cappé & Moulines, 2009). φ ∗n denotes the noisy estimate of the sufficient statistic for πn.
φ∗n = M1 M∗1 ∑ m∈Dk γ(δ (n)m ), (12a)
φ (k)n = (1−ηk)φ (k−1) n +ηkφ∗n , (12b)
πnewn = φ (k)n +αn−1
M1 +αn +βn−2 , n = 1,2, · · · ,N. (12c)
In terms of parameter w, the above practice is infeasible due to its non-exponential marginal distribution. Inspired by the stochastic gradient EM algorithms in Cappé & Moulines (2009), we perform the stochastic update in the original space. First, a local optima regression weight w(k) can be estimated via iterative optimization over the mini-batch Dk, using L-BFGS algorithm. Then we interpolate between the local optima and former estimations to form a global approximation.
w(k) = L-BFGS(L (w),g(w),Dk) (13a)
wnew = (1−ηk)wold +ηkw(k). (13b)
The convergence issues of the proposed online generalized EM algorithm are the analogues of the discussion given by Cappé & Moulines (2009) for their stochastic gradient EM Algorithms. The existence of such links is hardly surprising. In view of the discussions in Section 3 of Cappé & Moulines (2009), the online update rule (Eq.13b) could also be seen as a stochastic gradient descent formula, namely wnew = wold +ηk(w(k)−wold) .

heading: 4 EMPIRICAL ANALYSIS
text: In this section, we demonstrate the reliability of the proposed BDrank (Eq.5) with EEG signals from forty four participants.
Data Preprocessing: EEG preferences for each participant has been generated as follows: (1) the trials of each participant were randomly divides into two parts: 50% for training and 50% for test; (2) EEG preferences were constructed according to the pairwise comparisons between the RTs. To be specific, the type-1 preference propositions D1 were constructed with RT comparisons (RTm,1,RTm,2), which satisfies RTm,2 < min(RTm,2+τ1,τ2×RTm,2)< RTm,1; the type-2 preference propositions D2 were constructed with RT comparisons (RTm′,1,RTm′,2), which satisfies RTm′,2 < RTm′,1 < min(RTm′,2 + τ3,τ4×RTm′,2). It is notable that τ1 > τ3 > 0 and τ2 > τ4 > 1 control the difference in the RT comparisons simultaneously, we empirically set τ1 = 1;τ2 = 1.5;τ3 = 0.8;τ4 = 1.3 for all participants in our experiment setting. Considering the time delay among the channels in the time domain, Fourier transform (Welch, 1967) has been applied to EEG signals to transform timeseries into frequency domain. Further, to avoid overhead computation, EEG power within 0-30Hz has been selected, which is considered to be the most relevant to the RTs (Huang et al., 2015).
Baselines and Metrics: We compared BDrank with one regression model Support Vector Regression (SVR) (Smola & Schölkopf, 2004) and one ordinal classification model Logistics Ordinal Regression (LOR) (Harrell, 2001). First, we aggregate the predictions from different channels using a simple voting scheme, namely ρ̂m = sign ( ∑Nn=1 ρ (n) m [ 1(πn > κ)− 1(πn < 1− κ) ]) . ρ(n)m denotes the predicted state (1 means win and -1 means loss) for (RT 1m ,RT 2m) by the n-th channel, using the brain dynamics preference (x1n,m,x2n,m). ρ̂m is the final estimated order for (RT 1m ,RT 2m) by aggregating the predictions ρ(n)m over all channels. 1() is an indicator that returns one if the argument is valid and returns zero otherwise.
Then, we introduce two metrics to measure the performance of BDrank model from different perspectives. First, we adapted the Wilcoxon-Mann-Whitney statistics (Yan et al., 2003) to evaluate the accuracy (in %, higher is better) over all pairs, namely 1M1 ∑ M1 m=1 1(ρm = ρ̂m). Further, we investigate the reliability of BDrank in terms of preserving the global ordering w.r.t RTs. Note that a totally ordered set could be equally represented by a fully directed graph, where the fully directed graph can be further encoded by its degree sequence. Therefore, we first collected the indegree sequences3 (Becirovic, 2017) of the constructed directed graph using the predicted RTs and then measured the discrepancy between the predicted indegree sequences and ground truth using the root-mean-squared errors (smaller is better). See Supplementary for the detailed description.
Note that we only trust the predictions from informative channels with reliability πn > κ or πn < 1−κ . κ is set to 0.85 for all participants in our experiment. In terms of SVR and LOR, consider the scale difference between the EEG signals between different channels, we train a SVR/LOR model for each channel and aggregate the results from different channels to calculate the final predictions using the majority voting scheme. Since there is no mechanism for SVR and LOR to evaluate the channel state, we trust all the channels by default. Further, we only calculate the two metrics on the type-1 preference propositions D1, since the state of the type-2 preference propositions D2 is hard to evaluate when do prediction.
Parameter Initialization: In terms of the weight w, we randomly initialized w in [−1,1]. In terms of the channel reliability πn, we assumed all channels were relevant to the task beforehand, and randomly initialized the channel reliability πn,∀n = 1,2, · · · ,N in [0.5,1]. The L-BFGS implementation was downloaded from Granzow, where we used the default parameters. In terms of the hyperparameters (µ,Σ) for w, we adopted the standard Gaussian distribution, namely w ∼ N(0,1). In terms of the hyperparameters (αn,βn), as we intended to eliminate the effects of noisy channels, we adopted a strong non-informative prior for πn, namely αn = βn = 20, ∀n = 1,2, · · · ,N, according to Bishop (2006). In terms of the maximum iteration number, we set MaxIter = 30 in our experiment to ensure the algorithm converged for each participant. In terms of the minibatch size R, we set the sample ratio to 10% and 5% for type-1 and type-2 preference propositions, respectively. The parameters for LOR were optimized with L-BFGS with the default setting for all participants. In terms of SVR, we resorted to 5-fold cross validation to find the best parameters for each participant.

heading: 4.1 EMPIRICAL RESULTS OF BDRANK ON BRAIN DYNAMICS PREFERENCES
text: The accuracies of LOR and BDrank for training and test EEG preferences are presented in Table 1. For SVR, we followed standard procedure and optimization for each participant to predict RTs with training and test EEG signals. We adopted the Wilcoxon-Mann-Whitney statistics to measure the accuracy of the predicted RTs w.r.t. ground truth (shown in Table 1).
From Table 1, we observe that: (1) SVR achieves the highest accuracies (100% accuracy) on the training EEG data for almost all participants. However, test accuracies for most participants are inferior to other models, 21 participants for LOR and 35 participants for BDrank out of 44 participants. This observation is consistent with our previous discussion that regression-based models are easily overfitting to the training datasets, especially when extreme values (RTs in our problem) exist. (2) BDrank shows significant improvements over SVR and LOR in terms of test accuracies.
3We only consider the indegree sequence because the indegree and outdegree of a vertex can be uniquely determined when the overall degree of each vertex is fixed. The indegree of vertex vt can be calculated as deg−(vt) = ∑m∈N1(vt )[1(ρ̂m = 1) + 0.5× 1(ρ̂m = 0)] +∑m∈N2(vt )[1(ρ̂m = −1) + 0.5× 1(ρ̂m = 0)], where N1(vt),N2(vt) denote the index set of the pairwise comparisons with the RT of trial t (vertex vt ) appearing in the first and second position, respectively.
In particular, the number of participants with test accuracy above 75% are 17 and 25 for SVR and BDrank, respectively. This observation is consistent with our statement that classification-based models, served as a regularized alternative for the regression-based models, can effectively circumvent the overfitting caused by non-smooth RTs and preserve the ordering corresponding to RTs. (3) However, simple classification-based extensions (e.g. LOR) do not enjoy the benefits, which only achieves comparable test accuracies with SVR (outperforming SVR on only 21/44 participants). It is because they are not robust to the noisy data, while it is a common phenomenon when dealing with EEG signals. (4) Note that, compared to the fine-tuned SVR for each participant, fixed parameters were adopted for all participants in the experiment of LOR and BDrank. It further verifies the superior robustness of ordinal classification models compared to regression-based models.
To further investigate the reliability of BDrank in terms of preserving the global ordering corresponding to RTs, we first collected the indegree sequences of the constructed directed graph using the predicted RTs and then measured the indegree discrepancy between the calculated indegree sequences and the ground truth using the root-mean-squared error (RMSE) (shown in Table. 2).
From Table. 2, we observe that: (1) On the training datasets, SVR could exactly recover the indegree sequence (RMSE equalling to 0), which means SVR preserves the global ordering corresponding to RTs; while LOR and BDrank could not fully fit the training datasets. (2) Note that, on the test datasets, BDrank shows surprisingly good generalization performance and maximum preserves the global ordering corresponding to RTs.
To further demonstrate the superiority of our BDrank, we calculated the indegree sequences using the predicted RTs for participants P19, P38, P42, P43 with the most representative performance (Fig. 3). Similar observations can be found for other participants.
From Fig. 3, we observe that: (1) On the training datasets, the indegree sequences predicted by SVR completely overlaps with the ground truth; while the indegree sequence predicted by BDrank fluctuates locally but retains the overall trend. (2) On the test datasets, the indegree sequences predicted by BDrank still closely aligns with the ground truth with slight fluctuates; while the indegree sequences predicted by SVR fluctuates significantly and fails to maintain the trend with the ground truth, e.g. P19,P43. (3) It is interesting to note that the indegree sequences predicted by SVR usually fluctuates heavily for low indegree trials (denoting small RTs) and high indegree trials (denoting large
RTs). It means that SVR over-estimates the RTs with small values and under-estimates the RTs with large values. It is consistent with our claim that regression-based model is not suitable for the tasks with non-smooth response variable.

heading: 4.2 NOISY CHANNEL DETECTION
text: We also investigated the reliability of our BDrank from the perspective of noisy channel detection. According to our analysis, the parameter πn in the transition matrix Πn actually indicates the channel reliability. Hereafter, we leverage πn as the channel reliability indicator to detect noisy channels. Fig. 4 lists the noisy channels (marked in red) detected with 0.15≤ πn ≤ 0.85 , ∀n = 1,2, · · · ,N.
Fig. 4 shows that: (1) the noisy channels is universally exist among the EEG signals. There are total 42 among 44 participants with at least one noisy channels detected. For example, the 33-th channel is recognized as the noisy channel by BDrank for almost all participants. It is reasonable since the 33-th channel is the non-EEG channel, which is generally acknowledged as the non-relevant channel to any tasks; (2) For each participant, most channels are reliable, which ensures we can always find enough support to training our BDrank; and (3) The detected noisy channels varies from participant to participant, and do not possess the transitivity property between participants. Because the noise can arise due to (i) intrinsic non-EEG channel, e.g. the 33-th channel; (ii) channels for lateral mastoid references, e.g. the 23-th and 29-th channel (Chatrian et al., 1985); and (iii) improper experimentation or artifacts (Lin et al., 2018).

heading: 5 CONCLUSION
text: This work proposes a BDrank model to assess the state of mental fatigue. The efficacy of BDrank model was demonstrated using EEG data collected in sustained driving task from 44 participants. This model has been further combined with an online-generalized Expectation Maximization (OnlineGEM) algorithm to provide a continuous update in the model. BDrank model utilized a unique methodology with a regularized alternative, i.e. ordinal classification, to circumvent overfitting to the extreme values of RTs. It has been demonstrated that the overall performance of BDrank can be significantly improved with the introduction of the transition matrix, which enables the technique to evaluate the reliability of informative EEG channels while detecting noisy EEG channels. Empirical results show that BDrank combined with the OnlineGEM algorithm delivers significant improvement over SVR and LOR in terms of global ranking preservation.
In this work, the cooperation mechanism among channels is simplified as a weighted voting system, while different trials are viewed independently. We intend to further formulate it with more complex mechanisms, such as Markov decision process (MDP), to conduct learning and decision making simultaneously. Some previous (Chen et al., 2016; 2015) studied the decision making process among crowd (noisy) workers, which is promising to our setting to investigate the cooperation mechanism among noisy channels. Efforts are underway to apply this approaches in future work.

Review summaries:
The mental fatigue is an important factor in road accidents. Finding a direct mapping between EEG features and reaction time is difficult and error-prone, combining the noise measurement of EEG and individual variation of RT.  The authors introduce a measure called BDrank based on partial ordering instead of regression. Formulating the measure as a MAP problem, the authors propose a generalized EM algorithm for prediction. An online extension, relying on iterative L-BFGS optimization over mini-batches. 

The paper proposes an algorithm for mental fatigue monitoring, relating a subjects' EEG signals to their reaction time (RT) during a simulated driving task, as an ordinal regression problem. The authors argue that RTs could be heavily skewed and/or non-smooth, making traditional regression approaches unstable due to outlier values. They propose a brain dynamic ranking algorithm,  BDrank, using a generalized EM algorithm to estimate its parameters, and compare it to support vector regression and Logistic Ordinal Regression, where they show improved performance by accuracy and root mean squared error (RMSE) over a database of 44 subjects.