The STEW (Simultaneous Task EEG Workload) (Lim et al., 2018) has been selected for experimental purposes. The dataset consists of raw EEG data collected from 48 subjects across 14 channels in two experimental conditions. In one condition, the EEG data was recorded from subjects in the rest state while not performing any mental activity. In the second condition, a multitasking SIMKAP test was presented to subjects, and EEG data was recorded. In both cases, a sampling frequency of 128 Hz was used with 2.5 min of EEG recordings utilizing the Emotiv EPOC EEG headset. Every recording contains 19,200 data samples (128 samples x 150 s) across the following 14 channels: AF3, F7, F3, FC5, T7, P7, O1, O2, P8, T8, FC6, F4, F8, and AF4. Additionally, a subjective rating was collected after each task execution whereby users rated their experienced mental workload on the scale 1 to 9. The rating was a likert scale with 1 = “very, very low mental effort”; 2 = “very low mental effort”; 3 = “low mental effort”; 4 = “rather low mental effort”; 5 = “neither low nor high mental effort”; 6 = “rather high mental effort”; 7 = “high mental effort”; 8 = “very high mental effort” and 9 = “very, very high mental effort.” The rationale of using a perceived mental workload scores, was a form of subjective validation to verify whether a subject indeed experienced an increase in cognitive load load while performing the SIMKAP condition as compared to the resting condition.

Applying a denoising pipeline is an important step to pre-process the raw EEG data and to remove noise from it to facilitate subsequent analysis. In detail, this process follows the Makoto's pre-processing pipeline (Miyakoshi, 2018) including:

The key pre-processing step is the application of ICA which is utilized to separate the 14 EEG signal sources into independent components for each subject. Fourteen components are generated and used to automatically find and remove artifacts without human intervention using part of the methodology described in Nolan et al. (2010). In detail, the criteria for identifying bad components includes the computation of the z-scores of each component's spectral kurtosis, slope, Hurst exponent and gradient median. Spectral kurtosis is a parameter in the frequency domain indicating the component's impulsiveness variation with frequency. The slope of a component represents its mean slope of the power spectrum over two-time points. The Hurst exponent, also known as a long term memory in time series, tends to measure the tendency of a component to either regress to its mean or to catch up with an upward/downward trend. The gradient median is the median slope of the component's time course. All of the components exhorting values above and below ranges “z − score±3” can be considered as artifacts since they are outliers and significantly different from all the others. The value of ±3 for a z − score was adopted from Nolan et al. (2010) as part of automatic outlier detection taken from the FastR method. Finally, the inverse ICA has been executed to convert the remaining “good” components back in the original neural EEG signal.

A baseline of initial parietal and frontal electrodes was adopted following the electrode locations from the 10-20 international system to form different alpha and theta clusters for analysis and comparison purposes. These electrode locations were cross-referenced with locations, naming notation and electrode availability from the Emotiv EPOC EEG headset. The initial electrodes that are selected from the frontal and parietal locations are indicated as S1 and S2 in Figure 1. Due to the limited availability of electrodes from the Emotiv EPOC EEG headset (highlighted in green in Figure 1), three frontal and one parietal clusters were constructed. In detail, cluster combinations and electrodes, together with the channel aggregation approach used, is shown in Table 1.

Generating band ratios from EEG channels over time follows the methodology used in Borghini et al. (2014). The computation of the alpha-to-theta and theta-to-alpha ratios was done utilizing the average power spectral density (PSD) values from the alpha band from the cluster c − α and the average PSD values from the theta band from clusters c1 − θ, c2 − θ and c3 − θ as outlined in Table 1. Three alpha-to-theta and theta-to-alpha ratios are set to take different clusters from frontal electrodes and one cluster from parietal electrodes. The computation of the band ratios are given as follows:

where, c − α and c_x − θ are the respective alpha and theta clusters (from Table 1), with e an electrode in them, and x a cluster among those using the theta band (c1 − θ, c2 − θ, c3 − θ). The combination of the clusters in Table 1, jointly with their individual use, led to the formation of the following possible mental workload indexes (configurations):

where, at-1=c-αc1-θ, at-2=c-αc2-θ, at-3=c-αc3-θ, ta-1=c1-θc-α, ta-2=c2-θc-α and ta-3=c3-θc-α. In this study, a 1 s non-overlapping sliding widow technique is employed to segment long EEG data, and for each window, an index of mental workload can be calculated.

Extracting high-level features from MWL indexes is crucial since it allows the finding of distinguishing properties that otherwise would not be possible if a raw index alone is considered. The extraction of such high-level features from the indexes defined in the set Equation 3 is executed using the TSFEL (Time Series feature Extraction Library) (Barandas et al., 2020). The advantage of using TSFEL is that it offers a wide range of statistical properties that can be extracted from multiple domains, including those from frequency and temporal domains. It is useful for identifying peculiar aspects of a signal and its specific properties such as variability, slope or peak to peak, just to name a few. Classes of extracted features span from the most well-known such as statistical/spectral kurtosis, mean and median of a signal, to less frequently employed features such as human range energy ratio, the estimator of the cumulative distribution function (ECDF), variability and peak-to-peak. The idea behind considering a large number of initial features was to assess their individual importance, and subsequently retain only the most informative ones by adopting a systematic feature selection approach rather than selecting them subjectively from intuition. Feature reduction can also facilitate model training in terms of required computational time. The selection criteria were based on the “SelectKBest” feature selection algorithm that ranks the features with the largest ANOVA F-value between a feature vector and a class label. The reason for choosing such an approach is that it offers a better trade-off in terms of accuracy, stability and stopping criteria in comparison to other feature selection algorithms such as SelectPercentile or VarianceThreshold (Powell et al., 2019). Determining the threshold for an optimal number of features is an iterative process of supervised evaluation of model performance with variable numbers of features. Initially, a model with all features was built and its performance metric in terms of accuracy was observed and in subsequent steps, the number of features was reduced by half iteratively as long as the model performance increased. The following iterative step with number of features that would indicate a decrease in model performance served as a stopping criteria. Finally, a Pearson correlation was computed to assess the correlation between selected features in order to reduce multicollinearity among them. Reducing multicollinearity of features is an essential step for retaining the predictive power of each of them. Using highly correlated features very often hamper model training. Experiments conducted by Lieberman and Morris (2014) indicate a correlation threshold of ±0.5 for optimal model performance.

The modeling and training process aims at learning classification models capable of discriminating self-reported mental workload scores from subjects (target feature), given the features extracted and selected in the previous step D (independent features). The mental workload scores were selected rather than the type of condition (“Simkap” or “Rest”) because we wanted a sensible indicator of mental workload, not a task load condition. In other words, a self-reported indicator of mental workload can be considered a more reliable representation of the user experience than a class representing a certain task load condition. This argument is originated by the fact that, in both task load conditions, users can experience any level of cognitive load. For example, a novice user can experience high mental workload for an easy task load condition when compared to an expert user. Similarly, a skilled user can experience moderate mental workload even while in a resting condition because of significant mind wandering. In research from Charles and Nixon (2019), a distinction between objective elements of the work (taskload) from the subjective perception of mental workload is outlined. Both taskload and subjective perception of mental workload can be mediated by operator experience or time constraint factors. Therefore, it is intuitive that task load conditions are not equivalent to mental workload experiences. In fact, on one hand, the former are strictly defined prior task execution, and are static, meaning they are immutable during task execution. On the other hand, the latter are unknown prior task execution and can change depending on a number of factors, for example including user's prior knowledge, motivation, time of execution, fatigue, stress among the others. To stress further, research has clearly shown that even the same person can execute a task, designed with a specific, static load condition (pre-defined task demands) differently at various times of the day (Hancock et al., 1992).

Additionally, to facilitate subsequent interpretation, we treated model training as a binary classification problem, mainly to use more interpretable evaluation metrics such as precision, recall, accuracy and f1 score. Therefore, the target feature range of 1–9 of the self-reported mental workload scores was mapped into two levels of mental workload, the “suboptimal MWL” and “super optimal MWL.” The split was adopted based on the assumption of the parabolic relationship between experienced mental workload and performance as outlined in Longo and Rajendran (2021). This split was done by aggregating the scores from 1 to 4, representing some degree of low mental workload (effort), into the “suboptimal MWL,” and all the scores from 6 to 9, for all of those supporting some degree of high mental workload (effort), into the “super optimal MWL.” All those scores rated five were discarded because they represent the neutral experience of mental workload.

The learning techniques chosen for achieving this aim are Logistic regression (L-R), Support Vector Machines (SVM) and Decision Trees (DTR). Many research works have considered these three learning techniques for continuous and more prolonged EEG recordings (Berka et al., 2007; Hu and Min, 2018; Doma and Pirouz, 2020). Logistic regression and SVM, as error-based learning techniques, are suitable for binary classification tasks (as in this work). On the other hand, as an information-based technique, decision trees are suitable for distinguishing important features by calculating their information gains during model training.

Due to the fact that a small dataset of 48 subjects was selected, then repeated montecarlo sampling for model training and validation is set in the following order:

A general rule of thumb implies a minimum of 1/5th ratio for each feature in the data to increase model accuracy (Friedman, 1997). Given the low number of training instances in each of the target classes (“suboptimal MWL,” “super optimal MWL”), the “curse of the dimensionality” problem is anticipated (Verleysen and François, 2005). Therefore, a strategy for generating synthetic data is adopted, which is based on the generation of statistically similar synthetic data that mimics the original data. For this purpose, the Synthetic Data Quality Score based on metrics like Field Correlation Stability, Deep Structure Stability and Field Distribution Stability (Gretel.ai, 2022) is adopted.

The Field Correlation Stability is the correlation between every pair of independent features (fields) in the training data and then in the synthetic data. These values' absolute difference is computed and averaged across all independent features. The lower this average, the higher is the correlation stability of the synthetic data. Deep Structure Stability verifies the statistical integrity of the generated dataset by performing deep, multi-field analysis of distributions and correlations. This is done by executing Principal Component Analysis (PCA) on the original data, and comparing it against that on the synthetic data. A synthetic quality score is created by comparing the distributional distance between the principal components found in the two datasets. The closer the principal components, the higher the quality of the synthetic data. Field distribution stability measures how closely the field distribution in the synthetic data mimic that on the original data. The comparison of two distributions is done using the Jensen-Shannon (JS) distribution distance given as:

where H(O) and H(S) are the Shannon entropy values for original (O) and synthetic (S) data respectively and H(M) is the sum of selected weights for probability distributions (π) and dataset probabilities (P) given as M=∑i=12πiPi. The lower the distance score on average across all fields, the higher the Field Distribution Stability quality score and consequently the higher the quality of the synthetic data generated.

The Synthetic Data Quality Score represents an arithmetic mean between field correlation stability, deep structure stability and field distribution stability. In this sense, the Synthetic Data Quality Score can be viewed as a confidence score as to whether scientific conclusions drawn from the synthetic dataset would be indistinguishable as if they were to be used in the original data. Synthesizing new data is performed using synthetic generators offered from Gretel.ai¹. The training process for the combined (original + synthetic) uses the same montecarlo sampling with the same steps as with original data outlined above. Randomized 70% of subjects is selected both from the combined (original + synthetic) for the“suboptimal MWL” and the ‘super optimal MWL” classes (dependent feature) for model training, the remaining 30% of combined (original + synthetic) subjects was kept for model testing and performing 100 iteration through these randomized splits. During model training, the data was normalized using z − score normalization given as z = (x − μ)/σ, where μ is the mean of training samples and σ is the standard deviation. The rationale for using z − score normalization is that it tends to minimize the mean (μ = 0) and maximize the standard deviation (σ = 1) for the normalized value and makes it suitable reducing extremely peak values in data, by transforming it in such a way that it's no longer a massive outlier.

A set of evaluation metrics were employed to assess the ability of the selected models to generalize on unseen data by learning from the training data. These metrics can be used to measure and summarize the quality of the trained models when tested with previously unseen data. For a binary classification problem, such as in the case, the evaluation of the models is dependent on True Positives (tp) and True Negatives (tn) which denote the number of positive and negative instances that are correctly classified. It can be also conducted with the False Positives (fp) and False Negatives (fn) that denote the number of miss-classified negative and positive instances, respectively. According to this, several metrics are used to evaluate the performance of the trained models. The accuracy metric measures the ratio of correct predictions over the total number of evaluated instances. Accuracy is represented as, Accuracy=(tp + tn)/(tp + fp + tn + fn).Precision is used to measure the positive instances that are correctly predicted from the total predicted instances in a positive class, given as Precision=(tp)/(tp + fp). Recall measures the fraction of positive instances that are correctly classified, Recall=(tp)(tp + tn). F-Measure or f1-score is the harmonic mean between recall and precision values represented as, f1 − score = (2·Precision·Recall)/(Precision + Recall). The proposed evaluation metrics are essential to assess the robustness of the selected models built upon high-level features extracted from the MWL indexes toward the discrimination of self-reported perceptions of mental workload. While precision refers to the percentage of relevant instances, recall refers to the rate of total relevant instances correctly classified by the model. The best model minimizes the value of (fp) in precision and (tn) in recall, and both come at the cost of each other since we cannot minimize both of them in one metrics. f1-score represent a harmonic mean of precision and recall and takes into account both metrics. Consequently, to bring hypotheses H1 and H2 on provable grounds, the f1 − score metric is adopted too.

Artifact removal is performed on each EEG signal for each of the 48 subjects separately for the “Rest” and “Simkap” task load conditions. The average number of ICA components removed from the EEG data associated with each subject is 1.61 for the “Rest” and 1.46 for the “Simkap” condition. The number of removed artifacts is within limits of the adopted methodology described in Nolan et al. (2010). Figure 2 depicts the removal occurrence for a total of 14 components across all 48 users for “Rest” condition and “Simkap” condition. As it is possible to see from Figure 2, at most, one ICA component per subject that is significantly different from the other components (±3 standard deviations) exists. These components are removed by zero-ing them, and the EEG multi-channel data is subsequently reconstructed by applying inverse ICA. Since at least one bad component was identified and removed for most subjects, it is possible to reasonably claim that some artifact has been removed from the EEG signal, thus facilitating the subsequent computations of the alpha and theta bands.

All high-level features are collected from the statistical properties of the mental workload indexes in various domains, including the temporal and frequency domains. The initial number of collected features are 210, and the exhaustive list is provided in the Supplementary Material accompanying this article. The ANOVA F-Value is computed for each of these features, and those with the highest value are retained for subsequent model training.

Since the SelectKBest algorithm requires an initial number of features, an iterative approach of feature inclusion during model training and the performance of the models with those features is assessed through its accuracy. Iterative optimal feature selection is performed by employing data from the original dataset. Figure 3 illustrates the convergence on the optimal number of features in relation to model performance grouped by learning techniques (L-R, SVM, DTR). This resulted in a reduced number of features that are kept for model training by employing to the dataset enhanced with synthetic data. Figure 4 shows the Pearson correlation matrix for the “Rest” and the “Simkap” states for the alpha-to-theta ratios for the case of index at-1 (as designed in Section 3.3).

Noticeably, most of the features are in the correlation range between −0.5 and +0.5, which contributes to reduce multicollinearity and thus potentially being all relevant and with high prediction capability (Lieberman and Morris, 2014). Figure 4 is an illustration of the results associated to a single mental workload index (at-1). However, results associated to the other indexes are mostly consistent with these, as it is possible to examine in the Supplementary Figures S1–S9 accompanying this article.

After the feature selection process, training of the models was conducted with Montecarlo sampling using Logistic Regression (L-R), Support Vector Machines (SVM) and Decision Trees (DTR) as described in design subsection E). Model training suffered from the “curse of dimensionality” issue since it comprised 48 subjects across only the seven selected features. The number of training instances is low compared to the number of independent features retained for modeling purposes. This is followed by the peak phenomenon of feature inclusion, where the number of features and their cumulative discriminatory effect is essential for the average predictive power of a classifier, which is data-dependent (Zollanvari et al., 2019). The initial model evaluation with test data on the original dataset did not reveal the accuracy of more than 60% for the standalone mental workload indexes built upon the alpha and theta indexes alone (c1 − θ, c2 − θ, c3 − θ, c − α). An accuracy of 70% for the mental workload indexes built upon the alpha-to-theta and theta-to-alpha ratios was observed. An in-depth analysis of the learning curves associated to the classifiers indicated a model underfitting and an inability to generalize from test data. Moreover, analyzing the spectral entropy of the mental workload indexes revealed a small variance variation, as can be seen from the boxplots of Figure 5. Small data variance subsequently increases the bias influencing the model's ability to generalize. Thus, as expected, synthetic data generation was applied for training robust models.

The input for data synthesis was the initial dataset comprised of 48 subjects and 150 data points (2.5 min of EEG data per participant split into segments of 1 s) for each of the indexes designed in Equation 3. Two synthetic datasets are created, one for the “Rest” and one for the “Simkap” task load conditions in order to retain the original dataset's intrinsic properties. Table 2 illustrates the overall synthetic quality scores for the mental workload indexes set in Equation 3.

Findings suggest that the overall synthetic data score is always above 87% throughout all the mental workload indexes selected for the comparative analysis. The synthetic quality score was measured in the scale (1–20)%-Very Poor, (20–40)%-Poor, (40–60)%-Moderate, (60–80)%-Good and (80–100)%-Excellent. This suggests that the quality of synthesized data is excellent, and in line with similar studies (Hernandez-Matamoros et al., 2020). Consequently, the synthesized data was of an additional 180 synthesized subjects with 150 data points (2.5 min of EEG activity split into 150 segments of 1 s) for each mental workload index. The final combined dataset with original and synthesized data is now comprised of 228 subjects with 150 data points for each mental workload index, as defined in set Equation 3.

The training of the models with Logistic Regression and Support Vector Machines learning techniques utilized the linear optimizer since it offers speed and optimum convergence on minimizing a multivariate function by solving univariate optimization problems during repeated training of the model (Fan et al., 2008). In the case of model training with Decision Tree, a Gini index was used to measure the quality of split during the model build.

The classifiers performance is shown in Figure 6. The evaluation metrics are shown across all mental workload indexes and are presented in descending order. The best classification accuracy results are observed for those models built with Support Vector Machines (SVM) and Logistic regression (L-R). In order to acknowledge the best learning technique, a two-tailed t-test between the three learning techniques and the employed evaluation metrics was performed. The results indicated no statistically significant difference between Logistic Regression (L-R), Support Vector Machines (SVM) or Decision Trees (DTR). This indicates the validity of the training approach adopted from the design, which means that no matter the learning technique adopted, the results across all applied evaluation metrics are the same. Table 3 illustrates the p-value significance levels of the t-test between evaluation metrics for each learning technique used in the study. The t-test was conducted with a threshold confidence value of α = 0.05.

Further analysis of mental workload indexes of the alpha-to-theta ratios (at − 1, at − 2) indicates better performance than their respective individual counterparts used for computing those ratios (c1 − θ, c2 − θ and c − α). For example, in the case of all learning techniques (L-R, SVM and DTR), first two alpha-to theta ratio indexes (at − 1 and at − 2) show better performance than their individual counterparts (c1 − θ, c2 − θ and c − α).

In the case of the theta-to-alpha mental workload indexes, this is also seen in the first two indexes (ta − 1 and ta − 2). Figure 7 illustrates the performance of band ratios (alpha-to-theta and theta-to-alpha) across all evaluation metrics, given as density plots for the case of Support Vector Machines (SVM). The density plots for all other learning techniques are available on the Supplementary Figures S10–S11. Table 4 outlines the significance levels of a two-tailed t-test between the alpha-to-theta and theta-to-alpha ratio indexes with indexes used to construct those ratios. A comparison analysis of models average performance between original data and those enhanced with synthetic data is shown in Table 5. An analysis of the number of electrodes across alpha and theta bands as given in Table 1 outlined in the design Section 3 we can see a higher number of electrodes in indexes c1 − θ and c3 − θ in comparison to indexes c2 − θ and c − α. To asses the impact of the number of electrodes in overall performance of the models, a cross-plotting between indexes at − 1 vs. at − 2 and at − 3 vs.at − 2 as well ta − 1 vs. ta − 2 and ta − 3 vs.ta − 2 is analyzed. Figure 8 illustrates this cross density plot comparison of performance between the alpha-to-theta and theta-to-alpha ratio indexes. Furthermore, a two-tailed significance test between these band ratio indexes (at − 1 and at − 3 vs. at − 2 as well as ta − 1 and ta − 3 vs. ta − 2) reveals a statistically significant difference. Table 6 presents the p-value significance levels for confidence interval of α=0.05. The p-values are with Bonferroni correction applied, resulting in a significance level set at α = 0.005.