Development of an expert system for the classification of myalgic encephalomyelitis/chronic fatigue syndrome

Study design, data, and compliance with ethical standards

In the open-access metabolomics dataset used in this study, the cohort consisted of 19 controls and 32 patients, all female, who were between age and body mass index (BMI) frequency-matched groups. The blood samples were collected in Ethylene Diamine Tetra Acetic Acid (EDTA) tubes. Then, plasma was separated from the cells through centrifugation at 500 g for 30 min and stored at a temperature of −80 °C for further analysis. After that, global metabolomics was performed using four ultra-high-performance liquid chromatography/tandem accurate mass spectrometry (UHPLC/MS/MS). This automated service allows for the precise measurement of hundreds of metabolites in a wide range of categories (Germain et al., 2018). A total of 42,432 data points from 51 participants were collected using the Metabolon^® technology to identify and analyze 832 different metabolites. Amino acids (177), carbohydrates (26), cofactors and vitamins (28), energy (10), lipids (353), nucleotides (29), peptides (33), and xenobiotics (176) are among the kinds of metabolites that this study examines. Information about the metabolomics data set, which is available as Supplemental File 1, contains information about the 83 sub-pathways that can be further separated into the eight super-pathways that Metabolon^® has identified.

The sample size required for this study was estimated with MetSizeR based on the PPCA model and calculated by setting the false discovery rate to 0.05. As a result, a minimum sample size of 28 patients in total with 14 patients in each group was estimated. Although it was challenging to find ME/CFS patients and healthy controls who satisfied the study’s inclusion requirements, the sample size was more than that predicted by MetSizeR (Nyamundanda et al., 2013), a technique used to estimate sample size in metabolomics investigations. The Institutional Review Board for Non-Interventional Clinical Research at Inonu University gave ethical permission to this study (decision no. 2023/4512).

Methods

The schematic representation of the methodology used and proposed in the research is as in Fig. 1.

Feature selection: Feature selection helps reduce the dimensionality of a dataset, making it easier to interpret and faster to process. It can also improve the performance of machine learning models by removing irrelevant or redundant features that can lead to overfitting (Chandrashekar & Sahin, 2014; Guyon & Elisseeff, 2003).

Analysis of Variance (ANOVA) feature selection is a statistical method used to identify the most important features (or variables) that are most strongly associated with the target variable. The method involves calculating the F-value for each feature, which measures how much the variation in the target variable can be explained by the variation in the feature. The higher the F-value, the more significant the feature is thought to be in predicting the target variable (Nasiri & Alavi, 2022).

The F-value is calculated as the ratio of the mean square between groups (MSB) to the mean square within groups (MSW): (1)${\displaystyle F=\frac{M{S}_{between}}{M{S}_{within}}}$

where: (2)${\displaystyle MSB=\frac{\sum {\left({\overline{x}}_k-\overline{x}\right)}^2{n}_k}{d{f}_{between}}}$

n_k denotes the number of measurements in the kth class, ${\overline{x}}_k$ denotes the mean of samples for each class, and $\overline{x}$ denotes the mean of entire samples.

The mean square within groups is calculated as: (3)${\displaystyle MSW=\frac{\left(\sum \sum {\left(x_{nk}-{\overline{x}}_k\right)}^2\right)}{d{f}_{within}}}$

x_nk represents the nth value for the kth class.

The degrees of freedom are calculated as: (4)${\displaystyle d{f}_{between}=K-1}$ (5)${\displaystyle d{f}_{within}=N-K}$

K and N are symbols for the number of classes and the overall sample size, respectively (Johnson & Synovec, 2002).

Extreme gradient boosting (XGBoost): XGBoost is a popular machine learning algorithm that belongs to the gradient boosting family of algorithms. It was developed by Chen & Guestrin (2016) in 2016 and has since become widely used in the machine-learning community due to its efficiency, scalability, and high level of performance on structured data problems (Cao et al., 2022; Ghaheri et al., 2023; Homafar, Nasiri & Chelgani, 2022).

The basic equation for XGBoost can be written as follows: (6)${\displaystyle {\hat{y}}_i=\sum _{k=1}^K{f}_k\left(x_i\right)}$

where ${\hat{y}}_i$ isthe predicted value for the ith instance, K is the number of weak models, f_k(x_i) isthe output of the kth weak model on the ith instance, and x_i isthe feature vector for the ith instance (Farzipour, Elmi & Nasiri, 2023; Nasiri, Homafar & Chelgani, 2021).

The weak models in XGBoost are decision trees, where each tree predicts the residual error of the previous tree. This approach is called gradient boosting, as the subsequent models try to minimize the gradient of the loss function concerning the predicted values (Ayyadevara & Ayyadevara, 2018; Chelgani et al., 2023). The learning objective function for XGBoost is a regularized version of the loss function, which helps prevent overfitting and improves generalization (Liu et al., 2021; Maleki, Raahemi & Nasiri, 2023).

The objective function that XGBoost minimizes can be written as: (7)${\displaystyle Obj\left(\theta \right)=\sum _{i=1}^{n}l\left({\hat{y}}_i,y_i\right)+\sum _{k=1}^K\Omega \left(f_k\right)}$

where θ denotes the model parameters, n is the number of instances, l is the loss function that measures the difference between the predicted and true values, and Ω is the penalizing regularization function for complicated models (Maleki, Raahemi & Nasiri, 2023) and is computed as: (8)${\displaystyle \Omega \left(f_k\right)=\gamma T+\frac{1}{2}\lambda {∥\omega ∥}^2}$

where γ and λ are variables that control the penalty related to the quantity of leaves T and the weight of each leaf ω, respectively (Nasiri & Hasani, 2022; Nasiri, Homafar & Chelgani, 2021).

Random forest: Random forest is an ensemble learning algorithm first proposed by Leo Breiman in 2001 (Breiman, 2001). The algorithm is based on the idea of combining multiple decision trees to create a more robust and accurate model (Rodriguez-Galiano et al., 2012).

Each decision tree is trained on a random subset of the training data and a random subset of the features at each split (Gong et al., 2018). This introduces randomness into the training process, which helps reduce overfitting and improve the model’s generalization performance (Alam & Vuong, 2013). At prediction time, each tree in the forest independently predicts a class or a numerical value, and the final prediction is obtained by combining the individual predictions, for example, by taking the majority vote (classification) or the mean value (regression) (Breiman, 2001).

One of the critical advantages of random forest is its ability to estimate feature importance, which can be useful for understanding the underlying patterns in the data and selecting the most relevant features for the task (Li, Harner & Adjeroh, 2011). The feature importance is calculated by measuring the decrease in performance (e.g., accuracy or mean squared error) when a particular feature is randomly permuted, which provides an estimate of how much the feature contributes to the model’s predictive power.

Support vector classifier: A support vector classifier (SVMs) is a machine learning algorithm that is used for classification tasks. The SVMs are based on SVMs theory, which Vapnik (1999) introduced in the 1990s.

The idea behind the SVMs is to find the best possible hyperplane that separates two classes of data points in a high-dimensional space (Khairnar & Kinikar, 2013). The hyperplane is chosen to maximize the distance between the hyperplane and the closest data points from each class (known as support vectors) (Pradhan, 2012). This distance is known as the margin, and the SVM/SVMs is often referred to as a maximum-margin classifier (Amarappa & Sathyanarayana, 2014).

To train SVMs with a linear kernel, the algorithm first identifies the support vectors, which are the data points closest to the decision boundary (Alam et al., 2020). Then, the algorithm finds the optimal hyperplane that maximizes the distance between the support vectors of each class (Luta, Baldovino & Bugtai, 2018). During the testing phase, new data points are classified based on which side of the hyperplane they fall on (Amarappa & Sathyanarayana, 2014). One of the advantages of using SVMs with a linear kernel is that it can work well even in high-dimensional spaces, where the number of features is much greater than the number of observations. However, it may not perform well when the data is not linearly separable. In such cases, a non-linear kernel can be used instead (Ghosh, Dasgupta & Swetapadma, 2019).

SHapley Additive exPlainations (SHAP): SHAP is a method for explaining the output of any machine learning model. It was introduced by Lundberg et al. (2018) and it provides a way to calculate the contribution of each feature to the prediction made by the model (Štrumbelj & Kononenko, 2014).

The SHAP method is based on the concept of Shapley values, which is a well-known concept in the cooperative game theory (Ekanayake, Meddage & Rathnayake, 2022). The idea is to calculate the marginal contribution of each feature to the prediction by considering all possible subsets of features that could have been included in the model (Fatahi et al., 2023; Li, 2022).

The main equation for calculating SHAP values is: (9)${\displaystyle {\varphi }_i\left(x\right)=\sum _{S\subseteq N\setminus \left\{i\right\}}\frac{\left|S\right|!\left(\left|N\right|-\left|S\right|-1\right)!}{\left|N\right|!}\left[f_{S\cup \left\{i\right\}}-f_S\right]}$

where:

ϕ_i(x) isthe SHAP value for feature i of instance x

N is the set of all features

S is a subset of features, excluding i

$\left|S\right|$ isthe number of features in subset S

$\left|N\right|$ isthe total number of features

$f_{S\cup \left\{i\right\}}$ is the output of the model when features in S are present along with feature i

f_S is the output of the model when only features in S are present (Lundberg & Lee, 2017).

This equation calculates the contribution of feature i by considering all possible subsets of features that could have been included in the model. It calculates the difference between the model’s output when feature i is present and the output when feature i is absent, averaging over all possible combinations of features that include feature i. This results in a measure of the marginal contribution of feature i to the prediction (Bi et al., 2020).

In practice, SHAP values can be calculated using various methods, including tree-based algorithms and kernel-based methods. These values can then be used to explain the output of a model and to identify which features are most important for making predictions (Mangalathu, Hwang & Jeon, 2020).

Performance assessment: In the presented study, individuals are classified as either having or not having ME/CFS based on the extracted metabolic biomarkers. Predictions can therefore be divided into four groups: True Positive (TP) signifies correctly identifying individuals with ME/CFS, True Negative (TN) denotes the accurate identification of individuals without ME/CFS, False Positive (FP) involves misclassifying individuals without ME/CFS as having the condition, and False Negative (FN) involves misclassifying individuals with ME/CFS as not having the condition. The following metrics are used to evaluate the classification results: (10)${\displaystyle Precision=\frac{TP}{TP+FP}}$ (11)${\displaystyle Sensitivity=\frac{TP}{TP+FN}}$ (12)${\displaystyle Specificity=\frac{TN}{TN+FP}}$ (13)${\displaystyle F1-score=\frac{2}{\frac{1}{Sensitivity}+\frac{1}{precision}}}$ (14)${\displaystyle Accuracy=\frac{TP+TN}{TP+FN+FP+TN}}$

Another performance metric utilized was the Brier score, which measures the accuracy of the model’s probability predictions. (15)${\displaystyle BS=\frac{1}{N}\sum _{t=1}^N{\left(f_t-o_t\right)}^2}$

In which f_t denotes the probability that was forecast, o_t represents the actual outcome of the event, and N, and N is the number of forecasting instances.