Ensemble Machine Learning Model Trained on a New Synthesized Dataset Generalizes Well for Stress Prediction Using Wearable Devices

Based on the prior experiments and reviewed literature (Table 2), publicly available datasets including WESAD [6] and SWELL [17] were utilized in this study. Additionally, since these datasets all included Empatica E4 sensor biomarker data, the Toadstool [20], UBFC-Phys [21], Non-EEG Dataset for Assessment of Neurological Status (NEURO) [18], Wearable Exam Stress Dataset (EXAM) [22], AffectiveROAD [19] and Multimodal Sensor Dataset for Continuous Stress Detection of Nurses in a Hospital [23] public datasets, which also are collected using Empatica E4, were considered. Table 1 provides a summary of all datasets considered. The AffectiveROAD dataset was excluded from our experimentation due to its labeling protocol. In this dataset, subjects were self- and observer-scored on a scale that limited conversion to a binary stressed/non-stressed indicator for model training and evaluation purposes. Subjects were scored while driving in inner city and highway scenarios, and it is not clear which scenario under normal circumstances would be considered more or less stressful.

Additionally, the Multimodal Sensor Dataset for Continuous Stress Detection of Nurses in a Hospital [23] was not included in our experimentation for similar reasons. In this dataset, thirteen potential scenarios could be marked by the subject as stressful by using the Empatica E4 event marker button. We found a substantial number of marked sections that exceeded the expected duration of the event, with a further number of short events with no potential for a cooling down period to separate the perceived stressful event from a subject baseline (non-stressed).

The SWELL, WESAD, NEURO and UBFC-Phys datasets were specifically labeled for stress, either using self-report scoring (UBFC-Phys) or periodic (SWELL, NEURO, WESAD, EXAM), where stressful situations were applied during specific time periods of the experiment. The Toadstool dataset contained sensor biomarker data recorded during periods of game-play where stress may have been perceived as high, but was not labeled. However, it was included to validate the XGB, ANN and ensemble models developed as part of this study. The EXAM dataset was also included for some of our experiments. This dataset was recorded during single sessions of 10 students writing mid-term and final exams, with exam scores provided, but no specific stress-based scoring included.

As stress response is not an instantaneous physiological reaction, overlap may occur between the labeled periods within the datasets and the biomarker data, as shown in Figure 1. This overlap can introduce miss-labeled biomarker data into a training set when building machine learning models and will likely further produce onerous predictions in new datasets during inference. It is therefor important to remove this overlap by excluding a short timespan prior to the stressor phase, and likely a longer period when the stressor phase completes to allow for a cool-down period during model training. Additionally, any periods labeled not clearly defined as stressed or non-stressed (baseline) should likely be excluded. In this paper, the following exclusions were applied to the datasets utilized:

Of the biophysical and biochemical markers normally used to measure stress [12], in this study we opted to focus on Electrodermal Activity (EDA), and Heart Rate (HR). Previous works have included Accelerometer (ACC) bio-marker data as a feature, however considering that at least one of the datasets (Toadstool) involved substantial subject physical movement as part of their experimental setup, this signal was excluded. Temperature (TEMP) and Skin Temperature (ST) was further excluded due to being missing in two additional datasets (SWELL, UBFC-Phys).

While not all datasets included an HR biomarker, this can be algorithmically derived using Heart Rate Variation (HRV) or Blood Volume Pulse (BVP), and where required was generated using a Python script utilizing the BioSPPy library [26].

The statistical programming language R [27] version 4.0.4 was used for analysis in this study. The Empatica-R library [28] was additionally utilized to read the raw Empatica E4 wearable device data into an R data table for processing. Next, sensor signals were converted to the same sampling frequency, as the Empatica E4 samples the BVP sensor data at 64 Hz, the EDA sensor data at 4 Hz, and HR sensor data in spans of 10 seconds. Stress labels within the datasets were then standardized to binary indicators with non-stressed set to zero (0) and stressed set to one (1).

Initial data exploration was performed on only the NEURO, WESAD, SWELL and UBFC-Phys datasets, considering the WESAD and SWELL datasets were utilized in prior work as listed in Table 2, and all four contained at least the EDA and HR biomarker data for a substantial percentage of the observations available. Not all biomarker sensor data were available for all subjects within the SWELL and WESAD datasets, and in these cases the subjects were excluded. This resulted in a total of 9 subjects from the SWELL dataset being included (of 25) containing a total of 157,739 observations, and 14 of the 15 subjects from the WESAD dataset, containing a total of 26,385 observations.

Statistical summaries were produced for each dataset, as detailed in Table 3, indicating substantial class imbalance across the datasets, with the NEURO and UBFC-Phys datasets containing 70.1% and 65.8% labels marked as stressed, respectively. In contrast, WESAD and SWELL have labels marked as non-stressed in 63.6% and 71.4% of total observations. The Table also shows the mean and Standard Deviation (SD) of the biomarkers across the four datasets. Substantial variation is further noted across each dataset for both the EDA and HR biomarkers.

Histograms (Figure 2) of each biomarker across the datasets show the EDA biomarker being significantly right-skewed. It should be noted that, the SWELL dataset contains on average 4% or more observations than the NEURO, WESAD or UBFC-Phys datasets, while the UBFC-Phys dataset contains the largest number of individual subjects, but each recording segment being significantly shorter compared to the SWELL, NEURO and WESAD datasets.

Box plots highlight this variation further (Figure 3 and Figure 4), with significant outliers occurring for the HR biomarker (Figure 4) of the UBFC-Phys dataset.

In addition to box plots and statistical summaries, plots were generated for each of the four datasets to investigate correlation between the HR and EDA biomarkers and the labeled, binary stress metric (Figure 5). Correlation was consistently weak across all the datasets, with the lowest correlation observed in the UBFC-Phys dataset. Within the NEURO dataset, we noted higher correlation between the binary stress label and the HR biomarker, while correlation was low between the binary stress label and the EDA biomarker. This correlation was slightly higher for the SWELL and WESAD biomarkers. This difference in correlation across stress biomarker datasets can affect transferability of experimental results, and further limit the ability of machine learning models trained on any of these individual datasets to generalize well for new, unseen data.

One of the most crucial steps in building machine learning models is the choice of algorithm. A common approach for optimal algorithm selection involves conducting exhaustive grid-searches using a wide variety of algorithms and approaches, and determining the best potential algorithm based on the results achieved. By examining the reported results of previous experiments and their associated algorithms (summarized in Table 2), we find tree-based algorithms including Random Forest performing well compared to other algorithms, in line with findings reported by Mishra et al. [29], who noted that Random Forest is often the optimal model when using EDA, HR and RRI biomarker data as input features for stress-related machine learning models.

Gradient Boosting algorithms, compared to Random Forest, offer improved efficiency and predictive performance, at the cost of being prone to over-fitting. Random Forest builds a number of decision trees independently, where each decision tree is a simple predictor with all results aggregated into a single result. The Gradient Boosting algorithm, however, is an ensemble of weak predictors, usually decision trees. Additionally, in Random Forest, the results of the decision trees are aggregated at the end of the process, while Gradient Boosting instead aggregates the results of each decision tree along the way to calculate the final result. Popular Gradient Boosting algorithms include LightGBM [30], CatBoost [31] and XGBoost (XGB) [32].

For this study, XGB was selected as the primary machine learning algorithm due to its popularity, availability as an R package, and having been extensively utilized since its original publication in 2016 [32]. XGB provides a highly efficient implementation of the stochastic gradient boosting algorithm and access to a suite of model hyperparameters, designed to provide control over the model training process.

To test the generalizability of machine learning models built on datasets containing sensor biomarker data of a small number of subjects, two initial experiments were performed, as follows:

1. 1.

   Train Random Forest, SVM and XGB on SWELL, test on NEURO and WESAD with no additional feature-engineering, using EDA and HR biomarkers only.

2. 2.

   Train Random Forest, SVM and XGB on SWELL, test on NEURO and WESAD and generate additional features using statistical summaries.
   

The SWELL dataset contains the largest number of subjects and observations, nearly 6 times more than that of either WESAD or NEURO. Additionally, the SWELL dataset showed the lowest standard deviation for both HR and EDA biomarkers (Table 3), therefore, it was selected for training. However, we also conducted experiments where either of the smaller WESAD or NEURO datasets were used for training. As expected, these experiments yielded poor results due to the low number of subjects and observations leading to a lack of statistical significance. These results and their code are made available through the paper’s public code repository (see Supplement6.R).

The UBFC-Phys dataset was excluded as it contained 56 subjects, substantially larger than WESAD, NEURO or SWELL, which were used in the previous works that were reviewed. As the stress metric within the datasets was binary, logistic regression was selected as the XGB learning objective. Due to the metric imbalance within each dataset, class balancing was performed by using the scale_pos_weight parameter provided by XGB, which has the effect of weighing the balance of positive examples, relative to negative examples, when boosting decision trees. Random Forest and SVM models were included to create baseline and compare initial results to those reported in previous studies (Table 2).

Optimal hyperparameters for the XGB algorithm were identified using a grid search with 10-Fold cross-validation, and this was repeated for each experiment, to ensure the most suitable parameters were utilized prior to model training. The XGB algorithm was run for 5,000 rounds with early-stopping set to 3 rounds. This ensures that training stops when the evaluation metric fails to improve after 3 rounds of training in order to prevent over-fitting. For the initial two experiments, training was automatically stopped after 93, and 87 rounds, respectively. For the Random Forest model, the number of trees were set to 200 [33], while a radial kernel was selected for the SVM model with the cost parameter set to 5, which were found to be optimal based on a number of trial experiments which were run with varying kernels and cost parameters ranging from 1 to 10.

Prior to performing experiments where feature-engineering was applied, we investigated the dominance of biomarkers and found HR feature as a more important biomarker, relative to EDA, as shown in Figure 6.

Feature importance was determined using the built-in xgb.plot.importance function provided by the XGBoost R package. Here, importance is calculated based on overall gain, where gain implies the relative contribution of the corresponding feature to the model, calculated by taking each feature’s contribution for each tree in the model. A higher value of this metric when compared to another feature implies it is more important for generating a prediction.

Next, an additional 22 features were generated by calculating statistical summaries using a tumbling window approach at 25-second intervals. Intervals between 0.25 seconds and 60 seconds are common in the literature (see Table 2), and 25-second intervals were selected based on the highest level of correlation observed between the biomarkers and stress label, after a substantial amount of experimentation. Statistical summaries including the mean, median, max, min, standard deviation, variance, skewness and kurtosis of each biomarker were calculated as new features, with an additional feature for the covariance between the EDA and HR biomarker. The resulting variable importance after model training and evaluation is detailed in Figure 7 (top 10, in order of importance).

The HR biomarker and its derived features are again observed to be the dominant features, with EDA-based features serving a secondary role. These 10 features were utilized in all our experiments requiring feature engineering. The remaining 12 features were not utilized to prevent potential over-fitting to specific datasets, notably UBFC-Phys, with 56 subjects, and substantially shorter recording segments compared to the NEURO, WESAD and SWELL datasets.

In order to increase the statistical power of the underlying smaller datasets, all four smaller datasets, i.e. NEURO, WESAD, SWELL and UBFC-Phys were merged row-wise without randomization, to form the single larger dataset, StressData. This includes a total of 99 unique subjects resulting in more than 85% statistical power, constituting a total of 244,399 observations, as follows:

Two new experiments, numbered 3-4, were then performed on the StressData dataset, as follows:

Of the prior works reviewed, we found three instances of using Artificial Neural Networks (ANN). Li et al. [11] utilized a Deep Convolutional Network (DCNN) consisting of three hidden layers containing 32 and 16 units, respectively. This was followed by a Sigmoid activation layer, with max pooling applied at the end of each convolution. Predictive accuracy using the WESAD dataset was reported as 99.80%. Khan [15] tested Long-Short Term Memory Network (LTSM), Convolutional Neural Networks (CNN) and Bi-Directional Long-Short Term Memory networks, using both supervised and semi-supervised methods and found supervised models to perform better than semi-supervised models, with accuracy across predictions averaging higher than 90% on the SWELL dataset. Eskandar et al. [8] utilized a LTSM model fed with features extracted by using a CNN with a reported accuracy of 85% (F1 Score), noting difficulty on predicting the baseline relaxed state when using the NEURO dataset.

Based on recommendations from available literature [34] and in consideration of the different experimental protocols and stressed/non-stressed period durations contained within the datasets utilized, we opted for a standard feed-forward architecture as a baseline model, rather than an LTSM model, which is more suitable for time-series data. Our proposed network consists of three layers. The input layer contained 10 neurons to receive the 10 input features. The hidden layer included 5 neurons (half the input features), connected to a final linear layer.

Mean Squared Error (MSE) was selected as the loss function, a frequently used measure of the differences between values predicted by a model and the expected values:

where N is the number of observations, $\hat{y}$ is the predicted value and y is the expected value. Typically, binary and multi-class problems including XGBoost would utilize log-loss as a loss function instead:

where N is the number of observations, log is the natural log, y is the binary indicator (0 or 1), and $\hat{y}$ is the predicted observation. However, we found MSE to outperform log-loss by a factor of 10% for predictive accuracy, and therefor utilized MSE for the neural network model, while retaining the default log-loss for XGBoost. The precision, recall and F-score metrics were also calculated given the classification nature of the labeling of the SWELL, WESAD, NEURO and UBFC-Phys datasets.

The Rectified Linear Unit (ReLU) [35] activation function was applied after both the input and hidden layers. Training was performed on a dual-GPU system using the Keras and Tensorflow libraries for a maximum of 200 epochs using a batch size of 512. Early stopping was employed to prevent over-fitting if the validation loss stops improving for a duration of 5 or more epochs.

Ensembling [36] is a widely-used technique known to improve a decision system’s robustness and accuracy. The motivation for using ensemble models is to reduce the generalization error of predictions, as long as the base models are diverse and independent [37]. Algorithms such as Gradient Boosting and Random Forest utilize ensemble methods to combine the individual results of a large number of independent decision trees into a single prediction. This approach has been further extended to combine the predictions of unique, individual machine learning models via majority voting, weighted-scoring and other blending techniques to deliver more powerful and robust models.

To apply ensembling techniques to our own work, we blended the predictions from our XGB model with those from our ANN model using simple averaging. This resulted in four new experiments as follow:

Reviewing the initial results from experiment 7 when using XGB, ANN and their ensemble against the labeled metric for all 99 StressData subjects, we noted predictions for UBFC-Phys vastly outperformed those for WESAD, SWELL and NEURO across both XGB and ANN models. This was attributed to two potential factors; i) the substantially larger subject size of the UBFC-Phys dataset and ii) the substantially shorter recording segment size of UBFC-Phys samples.

We further noted a reduction in predictive accuracy in Experiment 8, when the WESAD dataset was excluded from StressData and used as a new, unseen validation dataset. This implies that LOSO validation alone may not be a good predictor of overall model generalization ability, likely due to an imbalance across both labeled metric and the significant variance in recording segment sizes across the datasets when merged into StressData.

To address this imbalance, we combined the EXAM dataset with StressData, resulting in a total of 129 subjects. Next, we split the stressed and non-stressed periodic segments of all 129 subjects into two groups. We noted the shortest segment across all subjects were from the UBFC-Phys dataset, at 201 seconds. We therefore further split each stressed and non-stressed subset across all subjects into smaller sets of 180 seconds each (3 minute intervals).

This resulted in a total of 3,758 samples of 3-minute segments, labeled as either stressed (2,962 samples) or non-stressed (796 samples). From these segments, random sampling can be performed to build a substantially larger balanced training dataset of any number and combination of segments. We named this new dataset SynthesizedStressData.

In addition to the 8 aforementioned experiments, we performed two final experiments:

Experiments 1 to 2 tested two approaches, the first excluding feature engineering, with the second including feature engineering. These were performed for measuring generalization of models built on a single dataset (SWELL) containing a small number of test subjects, but with large amount of observations. These models were then validated against two new, unseen datasets (WESAD, NEURO) containing a substantially smaller amount of observations.

The results indicate low predictive power with the best model achieving 68% accuracy using a 70/30 Train/Test split during training, with feature-engineering. More importantly, the results show that a substantially larger number of observations (from a small number of test subjects) within the training dataset did not assist in generalization.

Furthermore, all three models (Random Forest, SVM, XGBoost) reported a large number of false negatives (Type II errors) and failed to predict the stress state of the test subjects (Figure 8), with the SVM model performing slightly better compared to Random Forest and XGBoost.

Experiment 2 repeats Experiment 1 by utilizing feature-engineering. Accuracy, Precision, Recall and F1-scores improve slightly when predicting on the unseen WESAD dataset, but reduces for the NEURO dataset. Additionally, the resulting confusion matrix (Figure 9) shows all three models producing a larger number of false negatives (Type II errors), compared to no feature-engineering as with Experiment 1.

Once the datasets are merged to form the larger StressData dataset for training a machine learning model, predictive accuracy of 63% is achieved with no feature engineering (Experiment 3), increasing to 72.36% (Experiment 4) when using feature engineering, both using LOSO validation. Importantly, Type II errors have been substantially reduced with the introduction of a larger, merged dataset, with a slight improvement when using engineered features, as shown in Figure 10. As neither Random Forest nor SVM showed a marked increase in performance over XGBoost, and due to the efficient training speed of XGBoost, both Random Forest and SVM was excluded for Experiments 5 to 10.

Repeating Experiments 1 and 2 but using a weighted ensemble approach (XGB + ANN), predictive accuracy slightly improves to 51% (NEURO) and 70% (WESAD) for Experiment 5, with no feature engineering, then 58% (NEURO) and 68% (WESAD) for Experiment 6, with feature engineering applied, when training on only the SWELL dataset. Experiment 7 combines the benefits of a larger dataset (StressData) with feature engineering and a weighted ensemble approach, reporting predictive accuracy of 80.33% when using LOSO validation.

Figure 11 shows the numerical prediction plot for a randomly selected subject (S9) from the SWELL dataset (Experiment 7), with red detailing the stress period and blue indicating the ensemble predictions. The result shows good correlation with the period labeled as stressed. We note from 15,600 seconds onward, the subject is still showing an elevated level of stress, while the labeled recording period has been marked as non-stressed, likely not allowing for a cool down stage between stress and non-stress periods during biomarker recording.

Figure 12 shows the numerical prediction plot for a randomly selected subject (here W14) from the WESAD dataset (Experiment 7), with red detailing the stress period and blue indicating the ensemble predictions. This result shows the ensemble model reacting well when subjects are subjected to a period of stress during the experiment.

Comparing the results to those from Experiments 1 to 4, we note slightly better results for the ensemble (using LOSO validation) when both a larger, more varied training dataset is utilized in combination with feature engineering. Various combinations of weighting were tested and utilized in order to achieve the highest accuracy rates, with the XGB model generally outperforming the ANN in most experiments:

However, when training StressData with the WESAD dataset excluded (Experiment 8), predictive accuracy reduced to 59% when testing on WESAD as an unseen validation set. Various class balancing methods were applied to StressData including over, under and both-sampling across the labeled stress metric, with both-sampling producing the highest score of 59%, compared to over and under-sampling. This lower accuracy result implies that LOSO validation alone (Experiment 7) is not a good measure of model generalization, as the training process includes the full StressData dataset, apart from the single subject being left out per training and scoring iteration. This result further implies that training on disparate datasets merged without extensive class and feature balancing across both labeled metric as well as the recording segment size of biomarkers, may not result in a well-generalized model.

Next, in Experiment 9, the EXAM dataset was included into StressData, constituting a total of 129 unique subjects, forming SynthesizedStressData. Random sampling was then utilized for class balancing as described in section 3.8 to construct 200 training subjects consisting of two physiological states: 6 minutes of a non-stressed baseline, and 6 minutes of a stressed period, sampled from the pool of 3,758 3-minute StressData segments (2962 stressed and 796 non-stressed). Validation was performed using LOSO, achieving a mean accuracy of 89% on this instance of SynthesizedStressData.

Precision, Recall and F1 score metrics are higher compared to experiments 1 through 8. Recall represents the model’s ability to correctly predict the positives (stressed) out of actual positives (stressed), while Precision measures the number of predictions made that are actually positive (stressed), out of all positive predictions (stressed). In both cases, a higher score is desired. The F1 score is a good measure to use when seeking a balance between Precision and Recall. The prediction plot of subject X188 (from the SynthesizedStressData) is shown in Figure 13, with a near-perfect match when predicting for both baseline and stressed periods.

Adjusting ensemble weighting based on varying biomarker segment sizes (Experiment 7) during the prediction phase, is not an ideal approach. Therefore, this method was not utilized in Experiment 9, where a fixed weighting of 60% (XGB) and 40% (ANN) (similar to Experiment 5) was applied, given the strength of the XGB model over the ANN model in prior experiments.

Finally, for Experiment 10, the ensemble model trained on SynthesizedStressData was applied to the WESAD dataset to compare predictive power to those from Experiments 1, 2, 5, 6, and specifically Experiment 8. For this experiment, all WESAD data in StressData was removed and used as a validation set (similar to Experiment 8), resulting in a smaller total of 115 subjects used for training. Random sampling was again utilized to build a new instance of the SynthesizedStressData dataset as described in section 3.8, to construct 200 training subjects consisting of two physiological states: 12 minutes of a non-stressed baseline, and 12 minutes of a stressed period, sampled from the pool of 3,600 3-minute segments (2917 stressed and 683 non-stressed).

A 12-minute interval was used to closely emulate the experimental protocol used when building the WESAD dataset, specifically the stressed period (approx. 10 minutes of stress [6]). The random sampling, training and validation process was repeated for 10 iterations, with accuracy scores averaged across the 10 iterations. The results of these iterations are detailed in Table 4.

The ensemble method outperformed both the XGB and ANN models, with a mean accuracy rate of 85%, and a very low standard deviation (0.029) across all 10 trials. A fixed weighting of 45% (XGB) and 55% (ANN) was used, showing more stability, and closer to an ideal balanced 50%/50% weighting. This result, compared to those from Experiment 8, shows a well generalized machine learning model capable of predicting accurately on new, unseen data (WESAD).

Figure 14 shows a prediction plot of subject W15 from the WESAD dataset. We note a slight stress over-prediction during the baseline condition, with high correlation between prediction and label across the 24 minute experimental period.

Revisiting the findings from Experiment 1 through Experiment 4, we find a significant reduction in Type II errors as shown in the confusion matrix of all three models when plotted for test subjects W4 and W14 from the unseen WESAD dataset (Figure 15). While both XGBoost and Ensemble models perform equally well for test subject W14, we note the Ensemble model outperforming both XGBoost and Artificial Neural Network models for test subject W4.

To emulate a real-world scenario where a new, unseen dataset is provided with a potentially stressful experimental scenario, we trained and applied the full SynthesizedStressData (3758 samples, 129 subjects) ensemble model to the Toadstool dataset, using a balanced 50% (XGB) and 50% (ANN) weighting.

This dataset contains sensor biomarker data recorded during periods of game-play where levels of subject stress may have been perceived as high. No labeling is provided for this dataset to compare predictive accuracy. Figure 16 shows a plot of the ensemble predictions for subject T2. We note elevated levels of stress when game-play starts, with a period of stabilization, prior to an increase as game-play reaches the end stage. While no definitive conclusions can be made from this test, it does show correlation between predictions and perceived levels of stress in this particular experimental scenario.

The results from the ten aforementioned experiments are detailed in Table 5. When examining the difference between experiment 4 and 7, we noted higher precision, indicating the ensemble model is better capable of predicting the stress state when subjects are in fact, under stress, while the improved recall indicates a better measure of identifying true positives (stressed condition). For experiments 2 and 6, there is a reduction in precision and recall for the WESAD dataset (better for the NEURO dataset), indicating that the ANN model, when part of the ensemble, fared worse for this dataset, reducing the overall precision and recall compared to singular models used for experiment 2. These differences when predicting on new, unseen data steered us towards merging the datasets and ultimately generating a synthetic dataset to address these shortcomings of training on small datasets.

We also conducted an experiment where the large SWELL dataset was excluded from the SynthesizedStressData for training. As expected, this experiment yielded poor results due to the low number of subjects and observations leading to a lack of statistical significance. The result and the code for this supplementary experiment are available through the paper’s public code repository (see Supplement10.R).