Feature selection on movement imagery discrimination and attention detection

2.1 Experimental paradigm

As depicted in Fig. 1, each trial started with the presentation of a cross centered on the screen, informing the subject to be prepared. Three-seconds later, a cue was presented on top of the cross. An arrow that was unbalanced (Experiment 1—Fig. 2a and b) or balanced (Experiment 2—Fig. 2c and d) in the visual field pointed to either left or right on a computer screen to cue the subject on the left or right hand movement imagery tasks. The subject was instructed to perform the movement imagery during a 4-s period starting at the cue presentation. Then, both the cross and the arrow were removed from the screen indicating the end of the trial. The intertrial period was randomly jittered to be between 3 and 4.5 s long. Each experiment had 2 runs of 40 trials each.

2.2 Electrode settings

Data were acquired from 19 electrodes according to an extension of the standard 10–20 system (i.e., Fp1, Fp2, F7, F3, Fz, F4, F8, T7, C3, Cz, C4, T8, P7, P3, Pz, P4, P8, O1, and O2). All electrodes were referenced to linked earlobes.

Data were digitized at 256 Hz and passed through a fourth order 0.5–60 Hz band-pass filter. Each channel's raw EEG signal was epoched from the cue time point to 4 s after the cue. The presence of artifacts in the epochs was automatically detected through the commercial EEG software Brain Vision Analyzer, Brain Products GmbH, that checks for maximum allowed absolute value (50 μV) at any time point, or maximum allowed absolute potential difference (20 μV) between two consecutive time points. The epochs contaminated with artifacts (e.g., eye blinks, muscle artifacts) were excluded from further analyses. The EEG features were extracted from 50 to 75 epochs out of 80 that were not contaminated with artifacts.

Data were re-referenced in order to improve the spatial resolution for coherence estimates. Because of the limitation of each re-referencing method [27], the original data $\left(V_i^{\mathrm{REF}}\right)$ were referenced to a common-average and to a Laplacian filter $\left(V_i^{\mathrm{LAP}}\right)$, according to Eqs. 1 and 2. S_i contains the four electrodes surrounding the central electrode $V_i^{\mathrm{REF}}$. The Laplacian filter calculates a weight g_ij for each peripheral electrode which depends on the distances d_ij between the j peripheral electrodes in S_i and the central electrode i, according to Eq. 2.

2.3 Power ratios

Alpha (8–14 Hz) and beta (16–24 Hz) EEG frequency bands include rhythms that are reactive to movement imagery [32]. Alpha is an idling rhythm which is also termed “rolandic mu rhythm” when generated in a motor-related cortex area. Alpha amplitude decreases during the execution of, as well as with imagined, limb movement at motor-related cortical locations. Beta rhythm generally increases in amplitude at limb movement initiation and termination at motor-related cortical locations [32].

The epochs used to extract the EEG power ratios include the data from the whole imagery time period (0–4 s from the time the cue is presented). Five narrow frequency bands were defined: 8–12 Hz (low alpha band); 10–14 Hz (high alpha band); 16–20 Hz (low beta band); 18–22 Hz (mid-beta band); and 20–24 Hz (high beta band). The power in the frequency broad band 0.5–30 Hz was used to normalize the power contained in the narrow bands. The power spectral density (PSD) was calculated through a periodo-gram to provide high-frequency resolution (no significant differences were found between periodogram and Welch methods). Each band power was computed as the sum of all PSD components (in V²/Hz) in the corresponding frequency range. The PR feature matrices had 95 features (5 PR features × 19 electrodes).

2.4 Event-related potentials

Event-related potentials are slow, nonoscillatory EEG potential shifts in response to certain events (e.g., visual or auditory stimuli) [23]. The ERP in response to auditory or visual stimuli can be modulated in amplitude and latency by stimulus parameters such as intensity [23], spatial distribution [31], familiarity [11] and EEG phase [18], as well as attention [21].

Following PR extraction, the raw epoch 0–1 s (post-stimulus) was low-pass filtered at 4 Hz with an eighth-order Chebyshev type I filter and used for ERP extraction. This epoch demonstrated the best task discrimination in previous work [5]. The filtered 256 point time series was down-sampled to 10 points per second. The first eight time points of the down sampled time series represent the features to be extracted from each EEG channel. A feature matrix with 152 features (8 ERP features × 19 electrodes) was generated.

3.1 Cross-validation scheme

The cross-validation error used to evaluate the down-selection algorithms’ performance was calculated through a 10-fold double-loop cross-validation scheme. The optimal number of features to select (input parameter of the proposed algorithms) is optimized in the inner cross-validation loop. The performance of the classifier previously trained with the selected features, is validated in the outer loop. Therefore, 10 validation error values are calculated. However, the 10-fold cross-validation results are considerably affected by variability. For this reason, the whole procedure was repeated 10 times. The median of the 100 validation error values (10-fold × 10 repetitions) was defined as the estimate of the classification online performance using the selected features.

A new feature subset was calculated at every validation fold. The number of times that each feature is selected for validation is indicative of its relevance for discrimination. Likewise, channel discriminative ability is assessed by its frequency of selection. A channel is deemed selected when at least one of its particular features (either ERPs or PRs) is selected.

3.2 Linear discriminant classifier

A form of Fisher Discriminant Analysis that was shown to be robust for spatiotemporal EEG pattern discrimination was applied [28]. The canonical discrimination function z is the result of a linear transformation of the original data Y according to Eq. 3. The discrimination coefficients of the canonical discrimination function are denoted by the vector a:

The discrimination coefficients in a were calculated according to calculations described in the Appendix. The Y matrices were determined by the feature selection algorithms for each cross-validation fold.

The discrimination quality was assessed through three different measures: cross-validated error rate for group membership prediction of every p-variate sample; Wilks’ statistic and its normal theory confidence limits; and a bootstrapped confidence limit that is robust against deviations from normality in the data structure [10]. The multivariate data matrix Y is transformed to the vector z which has mean u and a normal p-variate distribution f(z). Prior probabilities π_j are calculated as the ratio of the samples within group j to the total number of samples n. According to the Bayesian theory, the probability that the data came from group j of two groups (in our case left-hand or right-hand movement imagery), given a vector z, is calculated through π_jz in Eq. 4:

The term π_jf_j(z) was approximated by exp[q(z)] where $q\left(\mathbf{z}\right)=u_j^{\mathrm{T}}\mathbf{z}-\frac{1}{2}{u}_j^{\mathrm{T}}{u}_j+\ln \left({\pi }_j\right)$ [10]. The mean of the canonical discrimination function for group j is u_j. The highest π_jz determined the predicted group membership for each sample. In order to robustly assess discrimination quality, the cross-validated prediction error rate was calculated. The normal theory Wilks’ statistic (W) analyses the eigenvalues of the linear transformation Ya above. Since W is chi-squared distributed, the discrimination significance was assessed by comparing W with the 95% confidence limit. The statistic W approaches the value zero as the group separation improves. As our data certainly deviate from normality, discrimination quality was also tested through a re-sampling technique (bootstrap) that permutes the data labeling to test whether the group assignment is meaningful. The confidence limit of a classification error was set to 95% of the 100 permutations that were tested. Further details on the discrimination quality measures can be found in [28].

3.3 Sequential forward selection algorithm

This method resulted in a wrapper-type algorithm since a feature is included in subset i if it leads to the highest group discrimination (lowest W_i) of the canonical function z_i among all remaining features. The four main steps of this algorithm are described below:

3.4 Across-group variance algorithm

This proposed filter-type method uses a special formulation of PCA [6] to select features while reducing data dimensionality. Initially, Y is decomposed through singular value decomposition (SVD) into three matrices: U_n×c (component orthogonal matrix; c is the number of principal components), S_c×c (singular value diagonal matrix), and V_p×c (eigenvector orthogonal matrix; p is the number of features). The eigenvalue vector λ of the feature covariance (Y^TY) is calculated as the diagonal of S². The principal components are linear projections of the features onto the orthogonal directions that best describe the data set variance. However, when data presents a group structure, the information provided by a component is more detailed than a variance value. Thus, the total covariance (ψ) can be decomposed into a sum of within (ψ_within) and between (ψ_between) group covariance parts. The pooled covariance matrix is used as an estimation of the within group covariance matrix (ψ_within) as in Eq. 5:

By using the Bessel's correction, the sample covariance matrix ψ_i is weighted by n_i – 1 (n_i is the number of samples belonging to the ith group) instead of n_i in order to correct the estimator bias (ψ_i has rank n_i – 1 at most). The variance information provided by a principal component in vector notation is deduced in Eq. 6:

where v_j is the jth eigenvector (a column of V_p×c matrix) and corresponds to eigenvalue λ_j. While ${\mathbf{v}}_j^{\mathrm{T}}{\Psi }_{\text{within}}{\mathbf{v}}_j$ is a function of the sample distances to their respective group mean, ${\mathbf{v}}_j^{\mathrm{T}}{\Psi }_{\text{between}}{\mathbf{v}}_j$ is a function of the distances between the respective group means. In the discrimination context, only the latter comprises useful variance information. Therefore, the distance between groups given by the ith component, normalized by its total variance, gives a relative measure to calculate the AGV as in Eq. 7:

The between group covariance matrix (ψ_between) is calculated fromψ – ψ_within.

Although the principal components are organized by decreasing order of total variance (eigenvalues λ_i), this order is optimized for orthogonality rather than discrimination between groups. Therefore, the components are ordered according to the across group variance (AGV) in order to take the data group structure into account.

The dimensionality reduction results from the truncation of the c principal components ranked by decreasing AGV. The truncation criterion is a cumulative sum percentage of the descending ordered AGV scores and was assigned threshold values (typically 60–90%) for variance truncation [16]. The optimal threshold value was found by cross-validation. If k components met the truncation criterion, the truncated component matrix U_n×k (k < c) is a lower dimensional representation of the original feature space and more suitable for group discrimination. The retained variance information was transformed back to the original feature space using a modified version of the spectral decomposition property as in Eq. 8. In order to determine the features which resemble the retained components with minimal information loss, an across-group covariance matrix (ψ_AGV) was calculated:

Note that AGV_i is used instead of λ_i in the spectral decomposition Eq. 8. Each diagonal value of ψ_AGV represents the variance of a particular feature accounted for by the k retained principal components and measures feature discriminability. A ranked list with all p features in descending order of discriminability was determined. Finally, the number (q) of features that determine the optimal feature matrix Y_n×q was optimized by cross-validation.

5.1 Algorithm comparison

The PRs were used for comparison of the proposed feature selection algorithms (SFS and AGV) with three other algorithms in common use (RFE, GA, and relief) to discriminate movement imagery responses. Table 1a (common-average reference) and b (Laplacian spatial filter) presents the PR classification errors for all tested feature selection algorithms and demonstrate that AGV algorithm performed better than the other algorithms. The AGV algorithm achieved cross-validation errors between 12.5% and 38.09%, selecting only 4 of 95 available PRs, on average. Although the SFS algorithm performed similarly to the RFE algorithm, the former only selected 8 (common-average) and 7 (Laplacian) PRs while RFE selected 50 (common-average) and 55 (Laplacian) PRs. The GA and relief algorithms ranked next in increasing error. The AGV algorithm achieved lower error values and selected fewer features than the SFS algorithm in most cases. However, while the SFS algorithm achieved discriminations that were significant through the Wilks’ statistic in every cross-validation fold, the AGV algorithm discriminations were less significant. Thus, SFS and AGV algorithms demonstrated particular strengths.

The feature selection of the proposed algorithms was evaluated in terms of EEG channels and PRs in specific frequency ranges. As illustrated in Fig. 3a and c, central channels C3 and C4 were frequently selected by both proposed algorithms to discriminate between left and right hand movement imageries. The parietal channels P3, Pz, and P4 were also relevant for the left versus right hand movement imagery discrimination when the AGV algorithm was applied. The frequency bands 8–12 and 10–14 Hz were selected by both algorithms for most of the validations, as illustrated in Fig. 3b and d.

5.2 Event-related potential discrimination

The proposed feature selection algorithms were also applied on the discrimination of ERPs in response to arrangements of the movement cues (left and right arrows). While five subjects agreed to participate in Experiment 1 (same data as in PR analysis), only three of them underwent the Experiment 2. The arrows in Fig. 2a and b cued subjects in Experiment 1. Although, as in Experiment 1, the cues pointed either to the left or to the right, they were balanced across the visual field in Experiment 2 (Fig. 2c and d). The cue changes were employed in order to detect possible perception related EEG responses [21, 31]. Figure 4 compares cross-validated classification errors from both experiments and demonstrates that Experiment 2 generally led to a classification error increase (except for Subject C with AGV algorithm). Classification performance degradation was observed for both re-referencing methods.

In Experiment 1, the introduced algorithms down-selected feature spaces from 152 to less than 13 ERPs with cross-validation errors between 14% and 33%. The selection of ERPs in response to left versus right hand movement cues is illustrated in Fig. 5 for Experiment 1 and in Fig. 6 for Experiment 2. Considering both re-referencing methods, the channels F3, C3, P7, P8, and O2 were most relevant for discrimination in Experiment 1, for Subject A. Additionally, parietal (e.g., P7 and P8) and occipital channels (e.g., O1 and O2) were often selected among subjects in Experiment 1. According to the results of AGV algorithm (Figs. 5a, c, 6a, and d), the selection frequency of channels P7, P8, O1, and O2 decreased in Experiment 2. Although the results of the SFS algorithm demonstrate a less accentuated decrease on the frequency selection of such channels between experiments, a decrease in selection specificity was observed in Experiment 2. While ERPs with latencies 200 and 400 ms were selected frequently in Experiment 1 (Fig. 5b, d), no latency was particularly selected in Experiment 2.