Exploiting Patterns for Handling Incomplete Coevolving EEG Time Series

Abstract

The electroencephalogram (EEG) time series is a measure of electrical activity received from multiple electrodes placed on the scalp of a human brain. It provides a direct measurement for characterizing the dynamic aspects of brain activities. These EEG signals are formed from a series of spatial and temporal data with multiple dimensions. Missing data could occur due to fault electrodes. These missing data can cause distortion, repudiation, and further, reduce the effectiveness of analyzing algorithms. Current methodologies for EEG analysis require a complete set of EEG data matrix as input. Therefore, an accurate and reliable imputation approach for missing values is necessary to avoid incomplete data sets for analyses and further improve the usage of performance techniques. This research proposes a new method to automatically recover random consecutive missing data from real world EEG data based on Linear Dynamical System. The proposed method aims to capture the optimal patterns based on two main characteristics in the coevolving EEG time series: namely, (i) dynamics via discovering temporal evolving behaviors, and (ii) correlations by identifying the relationships between multiple brain signals. From these exploits, the proposed method successfully identifies a few hidden variables and discovers their dynamics to impute missing values. The proposed method offers a robust and scalable approach with linear computation time over the size of sequences. A comparative study has been performed to assess the effectiveness of the proposed method against interpolation and missing values via Singular Value Decomposition (MSVD). The experimental simulations demonstrate that the proposed method provides better reconstruction performance up to 49% and 67% improvements over MSVD and interpolation approaches, respectively.

1. INTRODUCTION

Human brain is one of the most vital organs of humans, controlling the coordination of muscles and nerves. To communicate between human brains and computers, electroencephalogram (EEG) based on brain computer interface is provided.

EEG signals reflect electrical activities of a brain. EEG signals have various clinical and advanced scientific applications such as medicine, pharmacy, psychology, linguistics, and biology. Due to the usefulness of EEG signals, the study of brain electrical activities, through electroencephalogram (EEG) records, is one of the most important tools for the diagnosis of neurological diseases [1].

To have an accurate insight and improved understanding of the mechanisms causing widespread brain disorder, careful analyses of EEG records are necessary. Large amounts of EEG signal processing have been investigated recently for distinguishing epileptic seizures, emotions, and brain functions. For example, Guler explored the ability of desired and trained Elman recurrent neural networks, combined with the Lyapunov exponents on the EEG signals [1]. Polat proposed a hybrid system to detect epileptic seizure in EEG signals via two steps: feature extraction using Fourier transform, and decision making using the decision tree classifier [2]. Adeli et al. investigated discrete Daubechies and harmonic wavelets for the analysis of epileptic EEG records [3]. Besides, many investigators have developed different methods to better understand the dynamics of human brains through EEG analysis [4]-[7]. Since EEG signals are recorded from multiple electrodes placed on different locations along the scalp surface, the problem of missing values can be encountered frequently due to the disconnection of particular electrodes, lost signals, or the failure to report some of the measurements in time [8]. Unfortunately, most of the above approaches assume that the EEG signal is a complete data set as input. With the occurrence of missing values, they can cause inaccurate usage, distort results, or even degrade the performance of techniques.

Missing values are often simply discarded or ignored because they are deemed unimportant for analysis. However, if the presence of missing values is consecutive for a long period of time, the performance of methodologies on these dataset will not provide the accurate results. Hence, missing value estimation in EEG time series is necessary for the pre-processing in order to boost accuracy and reliable usage requirement.

On the other hand, the collected EEG dataset exhibits a large spatial and temporal relativity in their values. The size of the data is often a massively large matrix containing multiple dimensions. Since there are correlations among the variables of a large matrix expression of EEG signal over time, the necessity of variable selection is emphasized. It is not advisable to try to use the whole dataset for a process that would require only a small piece to be completed adequately. Therefore, it is better to mine meaningful variables from the large databases in order to save memory, processing power, and time as well.

To handle the above challenges successfully, we propose a new approach which satisfactorily fits on the available data when missing values are present. Our prime objective is to mine correlations and evolving behaviors of multiple electrodes by automatically identifying a few hidden variables, and then to exploit their dynamics for solving the problem of missing observation. Correlation implies that the observed dimensions of multiple electrodes are not dependent. Therefore, missing values can be deduced from others through hidden variables. Evolving behavior denotes that missing values can be estimated effectively based on neighbors’ observations of next time ticks and following their moving trends. To evaluate the effectiveness of the proposed method by considering accuracy, reliability, and complexity aspects, this paper demonstrates the performance of imputation for consecutive missing observations on two real different datasets of electroencephalogram (EEG) signals. The proposed method can be effective to capture a few hidden variables automatically as well as to compactly illustrate how to learn their dynamics for solving consecutive missing values. Moreover, its computational time scales up linearly with the duration of the sequences. We compare the performance of the proposed method to the interpolation and missing singular value decomposition methods known as MSVD.

The outline of the remainder of this paper is as follows; In section 2, we provide a summary of the available articles and research literatures related to time series with missing values. Section 3 describes the materials and proposed model setup. Section 4 illustrates detailed object experiments and experimental evaluations. Finally, in section 5, we summarize the whole paper and point out some possible future research directions.

 

2. RELATED WORK

In many real-world applications, time series have been given considerable attention within a variety of domains. They spread from network traffic data, currency exchange rates, sensor measurement, to biomedical and so on. In reality, the data set may contain missing observations since the process of data collection is not perfect due to poor record keeping, lost records, etc. Ignoring of missing data causes the loss of useful information of datasets. Therefore, numerous advanced filling mechanisms are proposed to overcome this problem in time series data.

One of the most straightforward procedures is to replace each missing variable with simple methods such as calculating appropriate mean values. Another alternating method for filling missing values is an interpolation method, which is related to the handling of missing elements using a curve fitting, known as linear interpolation and splines. These methods exploit the smoothness in a single time sequence. They can estimate the success for a short interval of missing values based on continuity of the sequence. Details of these approaches and its applicability can be found in [9]-[11]. However, these methods will either become invalid or face a big challenge when the observation gap of missing values is large. Furthermore, these approaches discard any relationship between the variables over time.

The well- known techniques related to dimension reduction and latent variables, namely Singular Value Decomposition (SVD) and Principal Component Analysis (PCA), could fill missing values through discovering correlations in multiple time sequences. They recover the possible missing values from estimation of other observations in one sequence of multivariate data. Many of the previous works has applied these methods for solving the presence of missing observations [12]-[14]. Kim at el. [15] proposes an incremental approach by updating the weights of PCA for every new data and the expectation-maximization to recover missing values. This model captured successfully correlations among different time sequences for imputing single missing value at various positions. However, this method ignores time-evolving patterns since its characteristics are not designed for tracking the ordering of the rows or the columns. Guha [16] identifies the correlations between multiple data streams by improving the discovery of correlations. They perform a dimensionality reduction with SVD periodically. However, this approach faces a challenge due to SVD re-computation and thus, it cannot easily handle the missing data. Papadimitriou [17] proposes SPIRIT model which can incrementally find correlations and hidden variables based on PCA and uses them to recover missing values. Nati and Jaakkola [18] provide a simple EM algorithm that factors a dataset to low- rank matrices and approximate the missing value problem. Olga Troyanskaya [19] proposes a method by using an expectation-maximization-like procedure. It performs SVD over all complete columns, regresses incomplete columns against SVD to estimate missing values, then re-factors and re-imputes the completed data until a fix point is reached. This is extremely slow and only works if only very few values are missing.

The recent intuitive mechanism for handling missing values by using a statistical model is known as the linear dynamical system, which can be used to estimate values for missing time points. Phong and Singh [20] demonstrate the efficiency of their method by exploring the linear dynamical system model for gene expression with missing values. Aristidou [21] uses Kalman Filter, which predicts the occluded marker position on human motion dataset. Dorfmuller and Kalman [22]-[23] show how to use the previous maker position and a skeletal model to estimate the missing marker locations using the extended Kalman Filter. However, these models become ineffective when the marker’s missing positions are held for a long period of time.

It can be observed from the literatures that there are a variety of techniques available for estimating missing values for time series data. However, only few literatures show the performance on different positions for the missing data, randomly in particular, for long-interval consecutive missing values. This research aims to propose a new approach to handle missing values automatically. Moreover, its computation time is significantly considered with the linear on the duration of the sequences.

 

3. PROPOSED METHOD

In this context, a brief introduction of the Linear Dynamic System (LDS), also known as Kalman Filter used to model multiple time series, is first presented in order to enhance how to find the hidden dynamics in multiple EEG time series signals. Then, a new proposed model is built to impute missing values in multiple coevolving EEG time series.

3.1 A Linear Dynamical System for EEG Time Series

A sequence of EEG is a multi-dimensional data since multiple electrodes are used to record the electrical activity along the scalp surface. A dynamical system can be modeled by a sequence of multi-dimensional EEG signals, denoted by Y = {y1, y2, y3..., yT} , where each vector yt denotes the data at each time ticks t = 1, 2, 3… T of dimensionality of M. This means that data from EEG time series can be presented by a matrix YM×T of the variables M and observed time ticks T. We consider EEG time series data are obtained from EEG signals in such a dynamical system. It builds a statistical model to represent the state of the hidden variables which are evolving to a linear transformation leading to the observed numerical time sequences. The model can learn the dynamics of the time series data [24], [25]. It captures the correlations among multiple electrodes by choosing a proper number of hidden variables. In particular, LDS for multi-dimensional EEG time sequence is modeled by the following equations;

where θ = {μ0,Q0, A,Q,C,R} is the set of parameters. μ0 is an initial state for hidden variables of the whole system. Vector yn and zn denote observed data sequences and hidden variables at time t, respectively. The transition dynamic matrix A relates to the transition of the state from the current time tick to the next time tick with noise {ωn}. Matrix C is the observation projection with the noise {εn} at each time t, meaning that the series of hidden variables zn are evolving over time ticks with linear transition matrix A. Moreover, the observed data sequences yn are generated from these series of hidden variables with a linear projection matrix C. All noise, ω0, ωi and εi (i = 1…T) are zero-mean normally distributed random variables with covariance matrices Q0, Q and R, respectively. In the model, only the observation of the system is presented. The state and all the noise variables are hidden. Overall, the definition and mathematical description of symbols used in the system is shown in Table 1.

Table 1.Definition and mathematical description of Symbol table

3.2 Proposed Model Setup in the presence of missing values

The problem of missing time ticks will be first formulated in EEG data by the proposed model system. In the time course experiment, consider a collection of M-dimensional EEG time series Y of a length T with lost measurements; the missing values of the observations are indicated by matrix W. The matrix W of missing observation is the same size as Y and is defined as below;

The time sequences are modeled based on LDS, as seen in the above equations (1) – (3), with an extra missing indication matrix W [26]. We utilize an expectation-maximization algorithm to impute missing positions through estimating their expectation of missing values, E[Ymiss|Yobs], conditioned on the observed values, where Ymiss and Yobs are the set of variables for the missing values and the set of the observed values in the sequence Y, respectively.

In order to handle the problem of missing values, the prime objective is to mine meaningful patterns via automatically identifying a few hidden variables so that their dynamics will be discovered to solve the problem of missing observations. This study focuses on exploiting the dynamic connectivity of the brain signals via two particular properties: namely, correlation and temporal continuity. To meet this problem, it needs to model the dynamics and hidden patterns of the observed time sequence by using sequences of hidden state variables Z. To model correlations, the model uses data sequences, which includes both observed and missing values, generated from a series of hidden variables via a M ×T linear projection matrix C at each time point, shown in Fig. 1 where H is the number of hidden variables. Therefore, if some of the values are missing, they are inferred from the hidden variables since its mapping automatically discovers the correlations among the observation dimensions.

Fig. 1.Architecture of Proposed Model Setup

On the other hand, to model the temporal continuity property, HxH linear transition mapping A is adopted, since the hidden variables are time dependent with the values determined from the previous time ticks. This means that matrix A is related to the transition of the hidden variables’ states over time, describing how the states move forward over time. Thus, the next time point only depends on the current time point. In this case, we first set an initial state for hidden variables at the begining time point with the set of parameters θ = {μ0,Q0, A,Q,C,R}. In the system, the joint distribution of Yobs, Ymiss, and Z is given by the following equation:

To achieve these goals above, the proposed model is given to find the optimal solution to maximize the expected log-likelihood of the observation sequence with respect to the model parameters θ = {μ0,Q0, A,Q,C,R}, the hidden variables and the missing observation E[Ymiss|Yobs].

In practicality, in order to achieve parameter estimation, it is necessary to find the maximum likelihood L(θ) = P(Yobs). Yet, it is known to be difficult to maximize the data likelihood in presence of missing values. Therefore, the expected log--likelihood of the observation sequence over the parameter θ is maximized by using Expectation-Maximization (EM) algorithm [27], which iteratively maximizes the expected complete log- -likelihood as in equation (6). To achieve the maximum likelihood estimation of the model parameters, the expectation–maximization (EM) method for learning LDS is utilized. The algorithm iterates between computing the conditional expectation of hidden variables through the forward-backward procedure in E-step and updating model parameters to maximize its likelihood in M-step for estimating missing values [25].

In summary, the proposed method is performed to achieve the best parameters θ = {μ0,Q0, A,Q,C,R} for the model. The method applied in this paper conducts three main steps: expectation, recovering missing values, and maximization. In more detail, the algorithm first guesses an initial set of model parameters in the expectation step. It uses Kalman filtering and Kalman smoothing to estimate the hidden variables based on observation and current parameters for each iteration. The general idea is to use a forwardbackward algorithm to compute the posterior expectations of hidden variables, E(zn|Y;θ), tick by tick based on the computation of the previous time tick. Given the data with missing values, the estimation finds the marginal distribution for hidden state variables after initializing missing values to be a random number with interpolation method. Both prior and conditional distributions in the model are Gaussian, thus, the posterior up to the current time tick is p(zn|y1,..., yT), which is also Gaussian given by:

We obtain the following forward-backward propagation equations. The values here are μn,Vn and Pn-1, given by:

The initial values are given by following equations:

Smoothing involves an initial forward recursion followed by a backward recursion. In the forward step, the values of the Kalman Filter equations are stored. In the backward step, these values are then used to initialize the Kalman smoother equations given by:

Smoothing involves an initial forward recursion followed by a backward recursion. In the forward step, the values of the Kalman Filter equations are stored. In the backward step, these values are then used to initialize the Kalman smoother equations given by:

The expectations are taken from the posterior marginal distribution p(zn|y1,..., yT) from the propagation of belief. Therefore, the expectations are obtained by the following equations:

In the recovering step, missing values are recovered by using Markov property in the graphical model from the estimation of hidden variables, (Fig. 1) with the following equations;

In the maximization step, the algorithm updates the parameter θnew by maximizing the expected log-likelihood using some sufficient statistics from the posterior distribution. To estimate the parameters, the expected log-like hood L(θ;Y) in equation (6), with respect to the components of θnew, is maximized. Taking the derivatives of equation (6) and making them be zero provide the following results:

On the whole, the proposed method for solving the problem of missing values in EEG time series can be summarized as follows;

 

4. EXPERIMENTAL EVALUATION

4.1 Dataset Acquisition

To demonstrate the effectiveness of the proposed method by considering accuracy, reliability, and the complexity aspects, we evaluate its performances in recovering consecutive missing values on two real different datasets of EEG signals. We compare our proposed method to the interpolation and MSVD methods. The study is examined on random ranges of multiple coevolving EEG time series with consecutive missing values using a variety of parameter settings over different real datasets.

The first dataset is the publicly available Epilepsy EEG database at http://epileptologiebon-n.de/cms/frontcontent.php?idcat=193&lang=3&changelang=3. Further details are found in the work of [28]. This dataset contains no missing values consisting of five classes denoted by A, B, C, D, and E. Classes A and B are composed of segments taken from the surface EEG recordings that were carried out on five healthy volunteers. Classes C, D, and E originated from the EEG archive of presurgical diagnosis. Class D was recorded from within the epileptogenic zone. Class C was from the hippocampal formation of the opposite hemisphere of the brain. While classes C and D contained only the activity measured during seizure free intervals, class E contained seizure activity. Each of the class contained 100 single channels with 4,097 time points. In this study, we utilize the two sets named A and E from the complete dataset.

The second dataset is taken from http://www.bbc-i.de/competition/ii/. This dataset was recorded from a normal subject during a no-feedback session. The subject sat in a normal chair with relaxed arms resting on the table and fingers in the standard typing position at the computer keyboard. This BCI (self-paced) dataset contained two classes. Class label ‘0’ is used for upcoming left hand movements and class labels ‘1’ is used for upcoming right hand movements. Each class contained 28 EEG channels in the following order: F3, F1, Fz, F2, F4, FC5, FC3, FC1, FCz, FC2, FC4, FC6, C5, C3, C1, Cz, C2, C4, C6, CP5, CP3, CP1, CPz, CP2, CP4, CP6, O1, O2 [29].

4.2 Performance Evaluation on Reconstruction Error

To conduct our research experiment, two aspects are considered to assess the effective performance of the proposed method against MSVD and the interpolation approach. The comparisons are carried out based on the estimation performance on different amounts of entries with missing values.On the other way, we evaluate this problem based on a fixed percent of entries with missing values with different average missing lengths. For each experiment setup, we created different positions of consecutive missing observation on random channels of BCI and Epilepsy EEG datasets. The experiments are repeated 10 times in order to avoid the random effect.We reported the average of the mean square error (MSE) in order to evaluate the quality of the proposed method. The MSE is calculated following the equation: where t denotes each time tick, is reconstructed data, and y mentions the input data.

Table 2.Reconstruction error over different rates of missing values 5%, 10% and 15%

For each dataset in each experiment, the Fukunaka’s principle rule is used as a tool to achieve a proper number h for the hidden dimension of the model by taking the Singular Value Decomposition (SVD) of the original data Y =U ⋅ S ⋅VT, where both U and V are orthonormal matrices, S is a diagonal matrix with singular values on the diagonal. To get the number of h, small singular values are typically set to zero. Therefore, we order the singular values, and then choose h at the one with the 98th percentile of the total sum of squared singular values.

The comparison is based on the difference among the three methods’ reconstruction errors at 5%, 10%, and 15% of missing values, as shown in Table 2. In all of the cases, both the proposed method and MSVD use the same number of hidden variables with 98% energy; the average length of the consecutive missing values is 35 time points. Table 2 demonstrates that in all different amounts of missing data over the range of 5%, 10%, and 15%, the reconstruction errors of the proposed method give the best results that have smaller errors than the interpolation and MSVD methods. Specifically, in the BCI (self-paced) dataset, the proposed method presents 0.00483 and 0.00733 average reconstruction errors, which are lower than those of the MSVD and interpolation, respectively. On this dataset, it shows approximately 77% and 84% improvement compared to that of MSVD and interpolation.

Similarly, in subject 2 of the BCI (self-paced) dataset, the performance of the proposed method also presents 75% and 85% improved reconstruction against that of MSVD and interpolation. In the epilepsy dataset, the perfromances obatin higher reconstruction errors compared to the BCI (self-paced) dataset since it is more complicated. However, the average MSE of missing imputaions also shows significant reconstruction with lower errors at 0.069 and 0.1495 over those of the MSVD and interpolation, and shows 49% and 67% improvement over that of the MSVD and interpolation on subject A. On epilepsy of subject E, Table 2 shows that the proposed method presents estimated values of 0.0702 and 0.1574, which performed 49% and 68% improvements in the average MSE compared to MSVD and interpolation.

Secondly, the comparision is based on the difference average missing length settings with fixed 10% entries missing values among the three methods’ reconstruction errors. Fig. 2 shows the efficiency of the proposed method based on a fixed 10% missing entries with different consecutive missing lengths, ranging from 10 to 100 consecutive missing observations for BCI (self-paced) dataset in Fig. 2(a), and from 10 to 80 consecutive missing values for Epilepsy dataset in Fig. 2(b), respectively. Again, the proposed method performs the best reconstruction among the three methods with MSE increasing slightly along with the increasing consecutive missing length.

Fig. 2.Average reconstruction errors versus different average missing lengths on BCI and Epilepsy.

To test the efficiency of the proposed method, Fig. 3 shows an intuition of qualitative reconstruction on channel 4 from subject 1 of the BCI dataset of the three approaches. We demonstrate different hidden variables, which correspond to 96%, 97%, 98%, and 99% of energy in the original dataset. In all of the cases, the best recovery is performed on the 98% energy with 100 consecutive missing time points. The proposed method performs 74% and 85% higher accuracy for reconstruction than MSVD and interpolation on the BCI (self-paced) dataset, respectively. In the figure, the dashed portion indicates original signals, the green line is the reconstructed signals by the MSVD method, and the red line depicts the reconstructed signals of our proposed method. It shows that our proposed method (red line) achieves the best reconstruction since it reaches very close to the original signals against MSVD and interpolation methods.

Fig. 3.Reconstructions vs. original signals

Since the proposed methodology captures correlations via discovering a few powerful patterns as well as learning dynamics from multiple coevolving EEG time sequences, our proposed method provides better results than the compared approaches.

4.3 Computation time performance

The computational complexity of the proposed method is shown in Fig. 4. It can be observed that wall clock times lie on almost a straight line over the duration of sequences. The computation time of the proposed method increases slowly with the input and the time duration T of the time sequences. In this experiment, we use 98% energy on all sets, and the learning step runs at the same number of iterations equal to 20; we also set 10% missing values of the original dataset for each run together with 30 average occlusion lengths.

Fig. 4.Computation time versus the duration

Although the proposed method presents good results in all cases of the experiments and linear computation time, it still gets a hard time because the parameters of the proposed method are learned by the EM procedure. It takes more time for iterations to converge since the computation time of forward-backward procedures need cubic time in the dimension of observations. Therefore, the time complexity of the proposed method will receive poor scaling with large dimensional data. We will improve the speedup of the proposed model in future study.

 

5. CONCLUSION

In this paper, we proposed a method to tackle the problem of consecutive missing values for multiple EEG time series, specifically on their reconstruction. The method solves consecutive missing values of real multi-dimensional EEG time series. It can also automatically discover a few useful patterns and learn their dynamics successfully in order to solve the problem with consecutive missing observations. In most cases, this method provides the best result compared to the alternative techniques such as interpolation and MSVD. On the other hand, performance results show almost a linear speedup as we increase the input of the data set.

In this paper, we have only examined the EEG time series data for missing value reconstruction. For further work, the proposed method will be examined with other EEG time series mining tasks such as compressing and forecasting. Furthermore, we will broaden our research to involve different types of data such as Electrocardiography (ECG) and similar multiple coevolving data in the time series model. Finally, we will continue our exploration on the theme of mining large coevolving sequences with the goal of developing fast algorithms.