Non spontaneous saccadic movements identification in clinical electrooculography using machine learning

Non spontaneous saccadic movements identification in clinical electrooculography using machine learning Roberto Becerra-Gar´ıa1, Rodolfo Garc´ıa-Berm´udez3, Gonzalo Joya-Caparr´os2, Abel Fern´andez-Higuera1, Camilo Vel´azquez-Rodr´ıguez1, Michel Vel´azquez-Mari˜no1, Franger Cuevas-Beltr´an1, Francisco Garcia-Lagos2, and Roberto Rodr´ıguez-Labrada4 1Universidad de Holgu´ın, Grupo de Procesamiento de Datos Biom´edicos (GPDB), Holgu´ın, Cuba idertator,afernandezh,cvelazquezr,[email protected], [email protected] 2Universidad de M´alaga, M´alaga, Espa˜na [email protected], [email protected] 3Universidad Laica Eloy Alfaro de Manab´ı, Facultad de Inform´atica, Manta, Ecuador [email protected] 4Centro para la Investigaci´on y Rehabilitaci´on de las Ataxias Hereditarias, Holgu´ın, Cuba [email protected] Abstract. In this paper we evaluate the use of the machine learning algorithms Support Vector Machines, K-Nearest Neighbors, CART decision trees and Naive Bayes to identify non spontaneous saccades in clinical electrooculography tests. Our approach tries to solve problems like the use of manually established thresholds present in classical methods like identification by velocity threshold (I-VT) or identification by dispersion threshold (I-DT). We propose a modification to an adaptive threshold estimation algorithm for detecting signal impulses without the need of any user input. Also, a set of features were selected to take advantage of intrinsic characteristics of clinical electrooculography tests. The models were evaluated with signals recorded to subjects affected by Spinocerebellar Ataxia type 2 (SCA2). Results obtained by the algorithm shows accuracies over 97%, recalls over 97% and precisions over 91% for the four models evaluated. Keywords: Saccade identification, clinical electrooculography, classification 1 Introduction The alteration of eye movements is one of the symptoms of many neurological diseases like Parkinsons syndrome, spinocerebellar ataxias or the Niemann-Pick disease [4]. Specifically in the Spinocerebellar Ataxia type 2 (SCA2) this alteration is an important clinical marker present in more than 90% of patients [29]. There are several kind of eye movements such as saccades, fixations and pursuits. Among them, saccades are critical to follow and evaluate subjects with SCA2. For instance, SCA2 patients have significantly slower saccades and with larger latencies than healthy subjects [29]. The analysis of this kind of movement is used very often in the researches conducted by medical community, hence its importance. A technique for measuring eye movements called electrooculography consists in capturing the electrical potential of the eyes to calculate its magnitude and direction. This technique is widely used in electrophysiologic tests [16]. The resulting signals of this recording process are named electrooculograms [6]. Exists several methods and algorithms for identifying saccades in electrooculograms, the vast majority of them based on kinetic thresholds [11, 14, 31, 26], using suppervised learning [28, 6], unsupervised learning [20] or other novel approachs [18, 22] like particle filters [8]. These methods were designed to work in a not constrained scheme having advantages in a lot of scenarios. They are usually evaluated against data from healthy subjects where the differences between saccadic and non saccadic movements are very evident. However, in electrooculography clinical tests these methods try to detect as many saccades as posible, not distinguishing which of them are spontaneous and which not. In a previous work [2], we proposed a method that identifies saccadic movements using a sample-to-sample approach. This method allows us to discriminate where a sample belong to a saccadic movement or not. Now, in this work we have the task to identify which of these movements are stimuli related using a feature-based approach. Here we set out to evaluate the use of machine learning algorithms taking into account the strengths of clinical tests of electrooculography to solve the proposed task. Our approach have to use only horizontal movement signals and stimulus signals, and do not require the use of thresholds or any other user input. To do so, a new set of features were selected for training the models taking into account characteristics of valid saccadic movements. To identify the ocurrence of saccadic movements we use an impulse detection method based on velocity thresholds. These thresholds are calculated adaptively with a modified version of the method proposed by Nystr¨om and Holmqvist [18]. Our algorithm uses a classification model for solving the presented task, so we evaluate four of them: Support Vector Machines [7], K-Nearest Neighbors [27], CART decision tree [5] and Naive Bayes [25]. The performance of the classification models were measured, obtaining very good results (>97% accuracy) in all of them. The rest of this paper is organized as follows: In section 2 we describe the designed experiments and available data. Section 3 is devoted to analize and comment the results. Finally, section 4 summarizes the main conclusions and future work lines. 2 Material and Methods To test the selected algorithms an experiment was designed. The first step was identify potential impulses and annotating them for building a labeled dataset. After, each classification method is evaluated with stratified k-fold cross validation. Finally, we compare the performance of the models using nonparametric statistical tests to select the fittest. Clinical tests of electrooculography are setup as follows. Subjects with their head fixed are seated in front to a monitor at a previously known distance. Then, they are commanded to follow a visual stimuli which appears and disappears from one side to the other in the monitor. Capturing eye movements in these conditions using electrooculography allows to researchers the identification of which saccades respond to stimulus and which ones are spontaneous. Also allows to calculate important features of these movements like latency, duration, amplitude, deviation and maximal velocity. The electrooculograms were recorded using the OtoScreen electronystamography device at a sampling rate of 200 Hz with a bandwith of 0.02 to 70 Hz. Records of 12 sick subjects with SCA2 were used for building a dataset with features extracted from signal impulses. Each one of the records have at least tests of 10◦, 20◦and 30◦of visual stimulation. Typically saccadic tests have at least one horizontal channel and one stimulus signal (Fig. 1). 0 5 10 15 20 25 20 15 10 5 0 5 10 15 20 Position ( ◦ ) (a) Horizontal channel 0 5 10 15 20 25 Time (seconds) 150 100 50 0 50 100 150 Velocity ( ◦ /seconds ) (b) Velocity profile Fig. 1. Typical electrooculography signal with 30◦stimulus angle of a subject which suffers SCA2. Red signals are the scaled stimuli signals. Blue signals are the horizontal channel (a) and its velocity profile (b) respectively. IPython notebooks [23] were used in conjunction with the Python language scientific facilities: NumPy [19], SciPy [12], Pandas [17], Matplotlib [10] and Scikit-Learn [21] for running the experiments. These notebooks provides a rich interactive environment wich eases the development of experiments and facilitates the interchange of information between the authors. The intention behind using Python powered technologies is that the resulting algorithm (including trained models) will be used at NSEog, a processing platform developed by the authors. 2.1 Signals preprocessing Before the identification of potentially saccadic impulses, two common tasks needs to be performed: denoising and differentiation. Noise removal is very important matter in order to eliminate non desired spectral components due equipment malfunction, poor analog filtering or biological artifacts. Differentiation allows to obtain the velocity profile used later by the algorithm. Median filter (Equation 1) has proven to be very robust in eliminating high frequency signal noise while preserving sharp edges. An study carried out by Juhola in 1990 demonstrated that this kind of filters is appropiate for eye movements signals [13]. To eliminate non desired noise present in the signals used in the experiment, we use a median filter with a window size of 9 samples (approximately 45 milliseconds) obtaining very good results. This is accomplished using the medfilt function of SciPy. yi=median{xj|j=i−k,...,j+k}(1) Due the discrete nature of these signals, numerical differentiation is employed to calculate the velocity profiles. According [3], Lanczos differentiators (Equation 2) with 11 points (N= 11) have good performance for signals with the same characteristics as the ones used in this experiment. f0(x∗)≈3 h m X k=1 kfk−f−k m(m+ 1)(2m+ 1), m =N−1 2(2) We implemented the rutine of a Lanczos 11 differentiator which have the following formula: f0(x∗)≈f1−f−1+ 2(f2−f−2) + 3(f3−f−3) + 4(f4−f−4) + 5(f5−f−5) 110h(3) 2.2 Impulses detection Saccadic movements are represented as impulses in a velocity graph as shown in Fig. 1b. Typically, this movements can be easily identified by its contrast in magnitude and shape with other movements such as fixations and microsaccades. However, for the same stimulus angle the range of values of true saccadic impulses vary from subject to subject. This situation is tied greatly on the degree of affectation present in the subject [24]. One of the critical parts of the algorithm is the identification of velocities impulses which can potentially be saccades. For that matter, a threshold is needed for knowing when the velocity has reached a certain value that can be considered as a saccade candidate. Due the inter-subject variability explained before, this threshold should not be fixed a priori. Also should be large enough to ignore in most cases other movements like microsaccades and fixations, and not too large to miss valid saccadic movements. 0 5 10 15 20 25 Time (seconds) 150 100 50 0 50 100 150 Velocity ( ◦ /seconds ) Fig. 2. Threshold estimated in a 30◦stimulus angle test of a subject with SCA2. Estimation of a good velocity threshold is not an easy task, depending also on the noise levels present in the signal. The method described in Algorithm 1 solved this problem by calculating the threshold in an adaptive way. This is possible by using a safety margin (σ) based on the standard deviation of the signal which get adapted to the noise level present in the signal. The safety margin (σ= 6) employed by [18] ignores too many valid saccadic movements in lower angle tests for subjects with SCA2. A value of σ= 3 seems to be adequate for most cases at the expense of the detection of more non valid impulses. Even when this have a penalty in runtime performance, the final accuracy of the method should not decrease significantly. Due the amplitude of this new impulses the classification model should avoid them. Algorithm 1: Modified version of Nystr¨om and Holmqvist [18] threshold estimation algorithm Input :velocity profile (Array of degree/seconds samples) Input :σ(Safety margin) Output: Threshold estimation begin velocities ←− Abs(velocity profile); last threshold ←− Mean(velocities)+σ*Std(velocities); current threshold ←− 0; while Abs(last threshold-current threshold)>1do selected samples ←− samples from velocities below last threshold; current threshold ←− last threshold; last threshold ←− =Mean(selected samples)+σ*Std(selected samples); return last threshold; Another problem with this method [18] is that requires an initial threshold. This adds a bit of subjectivity to the overall process because this value could vary from healthy subjects to affected ones. We introduce a little modification to the method which states that this initial threshold is calculated the same way as the thresholds inside the iterations in the Algorithm 1. The only difference is that our method uses all velocity samples to calculate the mean and the standard deviation. This approach seems to work very well for all analyzed tests. Once the threshold is calculated, the next task is found all the impulses above it. This is accomplished by finding samples grouped together that exceeds this threshold. The principle behind this algorithm is looping through the signal to find velocities above the given threshold. When we encounter with one of these points, we move to the left and to the right until the velocity is zero or cross it. This approach allows further refinement of the saccade start and ending points because the impulses usually get more samples beyond the real saccade limits. If the length of a detected impulse is not greater than 10 samples, then is discarded to avoid very small invalid movements. A typical output of this method is represented in Fig. 3. 2.3 Model evaluation Once identified the saccadic impulses candidates, we need to know if they are saccades and if they are related to the stimulus. For this reason, the strategy behind our approach uses human intuitive features to solve this task. To take advantage of the characteristics of the clinical tests, the following set of features was carefully selected: Angle: Integer denoting the amplitude of the stimulus which can take 3 values: 10, 20 or 30. Absolute Latency: Time between the start of the stimulus transition and the maximal velocity point of the impulse in milliseconds (ms). 0 5 10 15 20 25 Time (seconds) 20 15 10 5 0 5 10 15 20 Position ( ◦ ) Fig. 3. Identified impulses in the same signal used in Figure 2 Normalized Latency: Normalized version of absolute latency with values between 0 and 1. The value 0 means that the maximal velocity is in the start of the stimulus transition, and the value 1 means that the maximal velocity is at the end of the stimulus transition. Amplitude: Difference between the maximum value and minimum value in the impulse. Deviation: Relation between amplitude and the angle of the stimulus. Maximum Velocity: Maximum velocity achieved during the impulse in ◦/s. Maximum Acceleration: Maximum acceleration achieved during the impulse in ◦/s2. Maximum Jerk: Maximum jerk achieved during the impulse in ◦/s3. Direction: Take the value 1 if the impulse follows the direction of the stimulus or -1 in other case. End Relative Position: Values between 0 and 1. The value 0 means that the impulse ends in the left position of the stimulus and the value 1 means that ends at the right of the stimulus. Using the features previously selected, a dataset of signal impulses was created. For building this dataset, a human specialist aided by the NSEog classified the detected impulses in valid and non valid saccades (Figure 4). As results, 1797 valid saccades and 6809 not valid impulses were obtained, resulting in 8606 instances. Because we are using Python technologies, Scikit-Learn was selected as machine learning library, hence we are constrained to a restricted set of models im- Fig. 4. Impulses annotation with the NSEog platform. plemented in it. The main policy of model selection was family representation, meaning that we try to choose methods with different working principles. So we evaluate four different models: Support Vector Machines, K-Nearest Neighbors, CART decision trees and the Gaussian version of Naive Bayes. Support Vector Machines (SVMs) are a set of supervised learning methods very effective in high dimensional spaces [7]. There are also very versatile supporting a set of kernel functions. Scikit-Learn implements four kernel functions: linear hx, x0i, polynomial (γhx, x0i+r)d, rbf e−γ|x−x0|2and sigmoid tanh(γhx, x0i+r). Preliminary experiments shows that for the proposed task, the rbf kernel function have the best performance compared with the others. Further study are necessary to fine tune the parameter γof this kernel. K-Nearest Neighbors is a type of instance-based learning which can be used for supervised or unsupervised learning. Instead of creating a generalizing function, it stores all the data inside the models using different data structures like Ball Trees or KD Trees. The principle behind the algorithm is find a number of training samples nearest to the analized point and predict the label from it [27]. To train our model we tried several numbers of neighbors starting from 2, giving the best results when this value is equal to 3. The data structure used is determined automatically by the Scikit-Learn implementation using optimization techniques. Decision trees are nonparametric supervised learning techniques. This algorithm requires little preprocessing and its runtime performance is enough to handle real time tasks. This method split the data trying to infere decision rules which can be used to clasify instances. Scikit-Learn uses an optimized version of the CART tree whichs support classification and regression [5]. The implementation used here do not require any parameter by default. Naive Bayes classifiers are supervised methods based on Bayes theorem which assumes independence between every pair of features [25]. We used a gaussian version of this classifier implemented in Scikit-Learn. Like the decision trees, the default implementation of this statistical classifier do not require any parameters. As validation scheme we use an stratified 10-fold cross validation to evaluate internally the models. The metrics employed to measure the performance were accuracy (Equation 4), recall (Equation 5) and precision (Equation 6) [30]. The accuracy give us a general quality measure of the performance of the models, while the recall and the precision allows to know how well the model predict or miss predict valid saccadic movements. In the following equations, TP (true positives), TN (true negatives), FP (false positives) and FN (false negatives) are the items from the confusion matrix used to compute involved metrics. Accuracy =TP +TN TP +FP +T N +F N (4) Recall =TP TP +FN (5) Precision =TP TP +FP (6)