A noise-immune reinforcement learning method for early diagnosis of neuropsychiatric systemic lupus erythematosus

1.   Introduction

Systemic lupus erythematosus (SLE) is an autoimmune disease that manifests remarkable clinical heterogeneity. This disease affects multiple organs with variable severity ^[1,2,3], including joints, lungs, heart, skin, and blood cells ^[4,5,6]. Furthermore, it may cause serious organ complications, such as lupus nephritis (LN) ^[7,8]. Neuropsychiatric systemic lupus erythematosus (NPSLE), as the most unknown manifestation of SLE, has severe effects on the central nervous system (CNS), peripheral nervous system (PNS), and autonomous nervous system (ANS) ^[9,10]. The nomenclature of NPSLE published by the American College of Rheumatology (ACR) includes 12 CNS and 7 PNS manifestations, which is widely used in clinical practice and facilitates research ^[11].

Although recent studies have made great progress in understanding the pathogenesis of NPSLE ^[12,13], the early diagnosis of NPSLE is still challenging due to the lack of gold standards. To date, researchers usually diagnosis NPSLE through clinical manifestations and a set of tests to make assessments, including brain magnetic resonance imaging (MRI) ^[2,14,15]. Due to MRI being widely used, it is the choice of clinical evaluation of the work-up of NPSLE, which can reveal lesions in the brain tissue of the patients ^[16]. However, almost half of the clinically active NPSLE patients with the absence of MRI abnormalities because of the lack of remarkable imaging manifestations, which shows the limitation of MRI in the NPSLE diagnosis ^[17].

Magnetic resonance spectroscopy (MRS) is an alternative method for detecting the concentration of metabolites non-invasively in NPSLE. Although previous work based on MRS has devoted massive efforts, there is still a long way to improve the accuracy of early diagnosis of NPSLE. Firstly, due to the data collected by MRS usually being sample-limited, it is challenging to build a robust model in this condition. Besides, the distribution of samples on these metabolic characteristics is high dimensional and nonlinear, affecting the model's accuracy significantly. Recently, deep learning have been widely used in many domains, such as traffic flow forecasting ^[18,19], medical image processing ^[20,21], and intelligent systems ^[22,23]. Note that the deep learning technology requires high quality and a large amount of high-quality data ^[24]. In this paper, we employ multi-voxel proton MRS (MVS) instead of the standard single-voxel proton MRS (SVS) for more accurate diagnosis, considered 13 metabolites in 9 brain regions concurrently, and combined a total of 117 metabolic features in total. It is easy to fall into the local minimum value, so we believe this method does not apply to our dataset.

Support vector machine (SVM) ^[25] is naturally good at learning nonlinear dependencies with limited samples. Therefore, SVM is an ideal classifier to avoid the problem above. Unlike deep learning, SVM is more suitable for small and mid-sized datasets because its hyperparameters are determined by supported vectors. In theory, SVM can achieve the optimum solution in the global range. This algorithm has been adopted in many fields, such as computer vision and machine learning. Zhu et al. ^[26] combined SVM and YOLOv3 to grade food with an accuracy of 98.5% achieved by SVM while the accuracy of YOLOv3 is 85.7%, and the overall accuracy is 96.4%. Cai et al. ^[27] introduced a SVM-based ensemble learning framework for basketball outcomes prediction. Tan et al. ^[28] proposed a hybrid GA-SVR to effectively forecast the short-term traffic flow in a large scale. However, the application of SVM has shown that the appropriate selection of the kernel parameter $\gamma$ and the penalty factor $C$ is the key to improving the learning ability and generalization ability of the SVM. In this work, we employ the genetic algorithm (GA) to optimize the hyperparameters of SVM to learn the nonlinear dependencies in features. Although there are many optimization methods on SVM, such as Kouziokas et al. ^[29] introduced an SVM Kernel, which combined the neural network weight vector of particle swarm optimization (PSO) ^[30] with the SVM kernel of Bayesian Optimization for time series problems. Faris et al. ^[31] proposed an approach based on a nature-inspired metaheuristic called multi-verse optimizer (MVO) for selecting optimal features and optimizing the hyperparameters of SVM simultaneously. These works have been combined with parameter optimization, but few letters are available to analyze the NPSLE diagnosis.

In recent years, reinforcement learning (RL) has been widely used in the healthcare domain. Petersen et al. ^[32] used deep reinforcement learning (DRL) to find an adaptive and personalized treatment of multi-cytokine mediation sepsis. Yu et al. ^[33] studied the integration of causal factors into the RL process and proposed a causal policy gradient algorithm to obtain better dynamic treatment regimes (DTRs) in HIV management. Maicas et al. ^[34] proposed a method based on DRL to detect the breast lesions automatically from dynamic contrast-enhanced magnetic resonance volumes with high accuracy and less detection time. Inspired by the powerful capacity of RL, we explore the possibility of empowering RL to the early diagnosis of NPSLE to further improve the model's classification performance. In this paper, the multi-agent reinforcement learning (MARL) is used to find the optimal feature subsets from the subset space. The central idea is to "assign" an agent to each feature, which controls the individual features and decides whether the feature is included in the final feature subset. Furthermore, feature selection is another effective technique for improving SVM classification. Some of the features are noisy or irrelevant, which can increase the computation quantity or even reduce the performance of the classifier ^[35]. Therefore, effective feature selection is necessary to enhance generalization ability, avoid overfitting, and provide a better explanation ^[36].

Existing reports mainly focus on optimizing the hyperparameters of SVM, selecting feature subset, or applying it to other applications, for example, hospitalization expense modeling ^[37], land-cover classification ^[38] and Polystyrene binding peptides (PSBPs) identification ^[39]. In the sample-limited and high-dimensional situation, few studies related to the early diagnosis of NPSLE in the machine learning field. The objective of this letter is to optimize the feature subset and SVM hyperparameters. On the other hand, we explore whether the proposed method can improve the classification performance of NPSLE early diagnosis.

The main contribution of this paper is summarized as follows:

● We employ the multi-agent reinforcement learning to find the optimal feature subset.

● We propose a genetic algorithm optimized support vector machine classification model for early diagnosis of NPSLE.

● We conduct extensive experiments to demonstrate the feasibility of our model. The experimental results show the superior performance of the proposed model to state-of-the-art models.

The rest organization of this paper is as follows. Section 2 describes the dataset used and the background material related to this work. Next, we introduce and present the proposed model and its details in Section 3. The outcomes of the experiments are elucidated in Section 4. The discussion and conclusion are presented in Section 5.

2.   Materials and methods

In this section, we first explain the details of the dataset. Then, we briefly introduce the methods employed in this paper, including support vector machines (SVM), genetic algorithm (GA) and reinforcement learning (RL).

2.1.   Materials

We analyzed proton magnetic resonance spectroscopy ($^{1}$H-MRS) data from twenty-three NPSLE patients and nineteen age-matched volunteers as healthy controls (HCs) who were admitted to the 2nd Affiliated Hospital of Shantou University Medical College and Shantou Central Hospital. All volunteers had no history of hypertension or diabetes and no head trauma, neurological, psychiatric, metabolic diseases, or other systemic diseases that may affect the nervous system. All patients met the classification standard of American College of Rheumatology (ACR) criteria ^[3] SLE in 1997 and the 1999 ACR definitions for NPSLE. Because of the ethical implications, routine drug therapy for NPSLE patients was not discontinued in the case group during the study period. The mean age of NPSLE patients is 31.71 while the one of HC is 29.5; there is no statistical difference in age between the two groups ($p = 0.765$). This study is reviewed and approved by the Ethics Committee of Shantou University Medical College. All subjects voluntarily participated in the study and signed the informed consent.

In order to evaluate the changes of metabolites in different brain regions of patients with NPSLE, a 3.0T MRI with MRS was performed in these subjects. Nine areas of metabolic data were collected, including right posterior cingulate gyrus (RPCG), left posterior cingulate gyrus (LPCG), right dorsal thalamus (RDT), left dorsal thalamus (LDT), right lentiform nucleus (RLN), left lentiform nucleus (LLN), right posterior paratrigonal white matter (RPWM), and left posterior paratrigonal white matter (LPWM), and right insula (RI) in all subjects. And obtained thirteen metabolic features from each brain region, include Creatine (Cr), Phosphocreatine (PCr), Cr + PCr, NAA, N-Acetylaspartylglutamate (NAAG), NAA + NAAG, NAA + NAAG / Cr + PCr ratio, myo-Inositol (mI), mI / Cr + PCr, Glycerophosphocholine (GPC / Cho) + Phosphocholine (Pch), Cho + Pch / Cr + PCr, Glutamate (Glu) + Glutamine (Gln), and Glu + Gln / Cr + PCr. Nine brain regions and thirteen metabolic characteristics were combined into 117 metabolic characteristics used as features in the dataset and set HC as $0$ and NPSLE as $1$ as the dataset's labels.

2.2.   Support vector machine

Support vector machine (SVM) is a powerful tool for binary and multi-class classification model developed by Vapnik and Cortes ^[25]. The learning strategy of SVM is to maximize the interval, which can be formalized as the regularized hinge loss function minimization problem. The basic idea of the SVM classifier is to find a hyperplane with the largest geometric interval between two groups that can correctly divide the data set. This is a small sample learning strategy with a solid theoretical foundation. The final decision function is only determined by a few support vectors. The computational complexity depends on the number of support vectors instead of the dimension of samples, which avoids the "curse of dimensionality" ^[40] in a sense.

In this paper, we present the dataset illustrated in Section 2.1 by $D = \left\{ (x_i, y_i) \right\}_{i = 1, \ldots, N}$, where $x_i \in R^M$, $R^M$ is the sample space, $y_i \in \{0, 1\}$ are the corresponding labels. $N$, and $M$ are the number of samples and features, respectively. For each groundtruth $y_i$, SVM model is to predict the corresponding result $\hat{y_i} = f(x_i)$ according to $x_i$. $f(x)$ is the classification function of SVM, which is defined as Eq (2.1).

where $\omega^T$ is the weight vector of nonlinear mapped metabolic characteristics, and $b$ is a constant bias term. The target of SVM is to find the optimal hyperplane that is explained in Eq (2.1). The objective function is shown in Eq (2.2).

where $\frac{1}{2} ||\omega||^2$ is the regularization term that minimizes the $L_2$ norm of the weight vector $\omega$, and $C\sum_{i = 1}^N \xi_i$ is the training loss. Sepcifically, $C$ is a penalty parameter of the training loss, $\xi_i = \max{(0, 1-y_i(\omega \cdot x_i + b))}$ is the slack variable.

The dataset we produced is high-dimensional (117 features) with limited samples (39 samples). SVM is good at handling small and medium-sized, nonlinear, and high-dimensional classification and regression problems and has succeeded in many practical applications. To this end, a reliable classification model based on SVM is naturally used to diagnose NPSLE.

In order to map our nonlinear dataset to high-dimensional space, we use kernel functions. In this letter, the radial basis kernel function (RBF) kernel is selected because it can classify data with high dimensions, unlike the linear kernel. Besides, the RBF kernel function has fewer parameters to define than the sigmoid kernel and the polynomial kernel. Specifically, it only has one hyperparameter $\gamma$ that affects the "spread" of the kernel ^[38].

2.3.   Genetic algorithm

The Genetic Algorithm (GA) is to find the global optimal solution in the search space of complex problems, which is a computational model simulating the natural selection and genetic mechanism of Darwin's theory of biological evolution ^[41]. The genetic algorithm is simple to implement without a complex problem-solving process. Figure 1 shows the process of the genetic algorithm.

Firstly, designing a scheme to encode the given information into the specific bit string encoding and setting the random initial starting point as the "parents". The encoding improves the convergence speed of the algorithm. Then, constructing the fitness function and obtaining each individual's fitness. Thirdly, creating the next generation through genetic operations, including natural selection, crossover, and mutation, as shown in Figure 2. After a certain amount of reproduction, GA can get acceptable individuals that satisfy the termination criteria. Due to the survival struggle, the best choice of high probability population will reach the overall optimal level ^[42].

2.4.   Reinforcement learning

The basic model of reinforcement learning (RL) is the interaction between individual and environment ^[43], as shown in Figure 3. The individual/agent is the part that can take a series of actions and expect to get a higher reward $r$. The other part related to this called the environment. The whole process is divided into different time steps. At each moment, the environment and the agent will have the corresponding interaction. The agent can take certain actions, which are imposed on the environment. After receiving an agent's action, the environment will feedback the current state of the agent's environment and the reward generated by the previous action.

The goal of RL is to get more rewards in the long term. There is a mapping between state and action, generally expressed by $\pi$; that is, $a = \pi(s)$, which means a state $s$ can correspond to an action $a$. The agent executes the action $a_t$ in the state $s_t$ and then transfers to the next state $s_{t+1}$ with the probability $T(s_t, a_t, s_{t+1})$. At this moment, it receives the reward $r_t$ from the environment with a discount coefficient $\gamma$. So the sum of reward from the beginning of $t$ time to the training end $T$ is $R_t = \sum_{t^{'} = t}^{T} {\gamma^{t^{'}-t}}{r_{t^{'}}}$. $\gamma \in [0, 1]$ is the discount coefficient, representing every moment's reward obtained in the future will be reduced. In order to estimate the cumulative reward $R_\pi$ obtained by strategy $\pi$ in the future, the value functions including state value function $V(s)$ and state-action pair value function $Q(s, a)$ are used to quantify the cumulative reward in reinforcement learning.

if the expected reward of a strategy is greater than or equal to the expected reward of other strategies, then it is the optimal strategy $\pi^*$. The optimal state value function $V^{*}(s)$ and the optimal state-action pair value function $Q^{*}(s, a)$ can be defined according to the optimal strategy $\pi^*$ as follows:

in the process of RL training, the optimal strategy $\pi^{*}$ is found by solving the value functions $V^{*}(s)$ and $Q^{*}(s, a)$.

3.   Proposed method

First, we design MARL to select the optimal feature subset and then use the selected features to train the classifier optimized by GA. In general, GA for parameter optimization is related to chromosome design, fitness function construction, and system processing. GA designs a suitable chromosome through binary encoding, then computes the fitness value and evaluates it by the system architecture.

3.1.   MARL for feature selection

We transfer the feature selection into multi-agent reinforcement learning by assigning an agent to each feature. Figure 4 illustrates the structure of feature selection based on MARL. The structure consists of agents, environment, actions, state, reward and reward assignment strategy.

● Agents. For the dataset described in Section 2.1 with $117$ features, we create $117$ feature agents that decide to discard or select the corresponding feature.

● Environment. We regard the selected features, that is, feature subset, as the environment.

● Actions. Action "0" means that the feature is disabled and will be excluded in the feature subset while action "1" means the feature is included in the feature subset. Note that we employ the $\epsilon$-greedy ^[44] to balance the exploration and exploitation. Specifically, the agent chooses the best action based on what it already knows at most of the time while it moves randomly with the probability $\epsilon$.

● State. We extract the mutual-statistic information as the state of the environment. Assuming that the dataset after feature selection is the matrix $F = (F_{ij})_{N \times n}$, $n$ is the length of the selected features, and $N$ is the number of samples in our dataset, here is 39. We first extract the statistical information of $F$, including mean, median, standard deviation, maximum, and minimum, to construct matrix $F^{s} = (F^{s}_{ij})_{5 \times n}$. Next, we extract five types of statistical representation of each line in $F^{s}$, obtaining the mutual-statistic information $F^{ms}$. The dimension of the $F^{ms}$ is $5 \times 5$.

● Reward. For our task, it is essential to boost the model's classification performance, here is accuracy $Acc$, and decrease the number of the feature subset $n$. To solve this issue, we define the overall reward function as follows:

where $b$ is the upper bound of the size of selected features, and $a$ is the joint action. The above global reward function balance the size of the selected subset and the performance of the model.

● Reward allocation. In this paper, the problem extends to multiple agents, where each agent acts simultaneously. We need to coordinate the cooperation and competition between agents ^[45]. The communication between features is planned by reward allocation. Specifically, if an agent $i$ does not change its action at the current iteration, it will be excluded when allocating reward. The reward $R$ is assigned equally to the agent that participates in the current iteration.

Q-learning is a prevalent reinforcement learning algorithm. At the time stamp $t$, an agent performs an action $a_t$ in the state $s_t$, then the state of the environment turns into $s_{t+1}$, the corresponding reward $r_t$ is sent to this agent. In Q-learning, Q-value $Q(s_t, a_t)$ is the maximum long-term reward after acting $a_t$ in the state $s_t$, which observes the following rules:

where $\gamma \in [0, 1]$ is the discount coefficient. The pseudocode of feature selection based on MARL is detailed in Algorithm 1. After that, we obtain the feature subset $F$ that be sent to GA-optimized SVM.

3.2.   Chromosome design

The first step of GA is individual gene designing and coding. In optimizing the combination of the two SVM hyperparameters, the kernel function parameter $C$ and the error penalty factor $\gamma$ are real numbers. We adopt binary coding to transfer them into bit strings, encoding the binary string of $l_1$ and $l_2$ bits. The gene string structure is shown in Figure 5. By combining the binary code of $l_1$ and $l_2$ bits, the gene string is obtained.

Because the data input to SVM should be decimal, we need to decode it into a phenotype according to Eq 3.3

where $p$ is the bit string, $min_p$ and $max_p$ are the minima and maximum values of the parameter respectively, $dv$ is the decimal value of the bit string and $m$ is the length of the bit string.

3.3.   Fitness function building

The goal of GA is to optimize the performance of SVM, which helps SVM get better generalization ability. The training samples are divided into two parts, training data, and testing data. Because our data samples are very scarce and difficult to collect, we use k-fold cross-validation to prevent over-fitting. We choose the mean score of cross-validation as the objective function to be optimized, using Eq 3.4

where $k = 10$. Mean squared error (MSE) can evaluate the change degree of the data. The smaller the MSE value is, the better the accuracy of the prediction model to describe the experimental data is, which can be calculated by Eq 3.5

where $y$ is the groundtruth, and $\hat{y}$ is the predicted label.

3.4.   Selection operator

To ensure the evolution to the direction of optimization, we use the wheel-roulette selection operator, which is based on an individual's fitness value. The higher fitness it has, the higher probability it can be selected.

3.5.   Crossover operator

The crossover operator is used to select individuals selected by the selection operator to reproduce the next generation. Two different chromosomes exchange genes at the exact location to create a new chromosome. The crossover operator includes two steps: random pairing and crossover. This paper uses a two-point crossover method; two crossover points are set in a single code string, and then partial gene exchange is performed. Take the second coding bit in the gene string as the first crossover point. Then, randomly generate a crossover point in the remaining binary code part, and exchange the corresponding gene fragments between the two crossover points. Figure 6 is an example describes the two-point crossover.

3.6.   Mutation operator

The mutation operator increases the ability of the genetic algorithm to find the global optimal solution. The mutation operator randomly changes the value of a position in the string with a certain probability. Specifically, converting the value of the binary-coded gene string from $0$ to $1$, or $1$ to $0$.

3.7.   System processing

The structure of the SVM classifier based on MARL and GA is presented in Figure 7, and the detailed steps are described below:

1) Assign a feature for each agent and perform action "save" or "discard" through executing the optimal policy to select the optimal feature subset;

2) Initialize the population to generate a number of individuals, including the SVM hyperparameters ($C$, $\gamma$);

3) The two hyperparameters ($C$, $\gamma$) are substituted into SVM to train classifier;

4) Compute the fitness of each individual. The average score of cross-validation is taken as the objective function to be optimized, and $C$ and $\gamma$ are taken as decision variables;

5) Judge whether the ending condition is met. If the termination condition is met, exit the cycle and finish the genetic optimization to get the optimized hyperparameter combination($C$, $\gamma$). Otherwise, go to step 6;

6) The individual with a high fitness will be chosen and maintained in the next generation;

7) After the selection, crossover, and mutation operator are executed to form a new generation of individuals, return to step 3 to continue.

Furthermore, to reduce the influence of large numerical features on the experiment, we scale the data before the first step, or the model may perform poorly. Here we use the $robustScaler$, which is defined as follows

where $v_i$ is a sample value, $median$ is the sample's median, and $IQR$ is the interquartile distance of samples. To overcome the disadvantage of limited samples, we use stratified sampling five times. Compared with K-Fold, it can ensure that the proportion of samples in the training set and test set is the same as that in the original data set.

4.   Experiment

We use MARL for feature selection to eliminate features that may be irrelevant or redundant, then we obtain 27 features, as shown in Table 1. Table 2 lists the initial parameters of the GA. Choosing these parameter values is based on the highest fitness value of the population.

Except for cross-validation score, sensitivity, specificity, and accuracy are used to evaluate the performance of these classifiers, which are defined in Table 3. We compare the performance of the SVM and other models, including Decision Tree, Bagging, AdaBoost, SVM and their combinations with PCA, SVM based on VAE, SVM based on GA, and the combination of SVM, MARL and GA. As shown in Table 4, the sensitivity, specificity, accuracy, and cross-validation scores of SVM are 73.9%, 75.0%, 74.4% and 0.754, respectively, while the corresponding metrics of our proposed method are 91.3%, 100%, 94.9% and 0.870, respectively, which are much higher than others and achieve the best results. In addition, Figure 8 presents the receiver operating characteristic (ROC) curves of different methods. In the case of SVM based on GA and MARL, the area under the curve (AUC) is 0.969, which is much higher than the AUC of other models, including other deep learning methods, such as the combination of VAE and SVM. The AUC of SVM and SVM based on PCA are only 0.854 and 0.867, respectively.

Figure 9 shows the changes of optimal individual's fitness and the average fitness of the whole population. At the 18th generation, the population stops evolution because the optimal individual reaches the judgment threshold that stops the algorithm. In the final evaluation phase, the fitness value of the optimal individual has risen to 0.87.

In addition, we compare genetic algorithm with grid search, a popular parameter selection method. Through GA, we get the optimal value of $C$, and $\gamma$ are 240.37, 0.0078, respectively. Therefore, we set the values of $C$ to be selected are $[240.30,240.37,240.40]$, and the ones of $\gamma$ are $[0.0075, 0.0078, 0.0080]$. We explore the optimal parameter combination through grid search, and the result displayed in Table 5.

By taking different values of $C$ and $\gamma$, we can explore the optimal combination of hyperparameters. However, this method can be time-consuming, especially when there are many hyperparameters that need to search. Besides, we set the value of $C$ and $\gamma$ above in the condition that we get the optimal hyperparameters through GA. If the range of values we set manually is unreasonable or the interval between different values is too large, it will be challenging to find the optimal value. By contrast, GA can search the optimal parameters automatically.

5.   Discussion and conclusion

In this paper, we propose a method that combines the feature selection based on MARL and parameter optimization based on GA to improve the early diagnosis of NPSLE. The data obtained by multivoxel proton magnetic resonance spectroscopy is high dimensional, sample-limited, and non-linear. It is hard to meet the model's requirements through traditional modeling methods or linear dimensionality reduction. SVM is an ideal classifier for our study, which has an advantage at dealing with small samples ^[46,47] and high-dimensional datasets. Even if a limited training set produces some bias, the model still has good generalization if we select appropriate parameters ^[48]. Besides, since the original data contains a lot of redundant information, it is essential to select the proper subset of features to avoid reducing the classifier's performance and increasing the computation quantity.

The proposed approach results were compared with Bagging, Decision Tree, AdaBoost and SVM without optimization to determine GA's ability of optimizing SVM hyperparameters. Experimental results show that the GA-SVM is significantly better than those classifiers just mentioned, achieving satisfactory sensitivity, specificity, and accuracy as 91.3%, 93.8% and 92.3%, respectively, while the SVM with default parameters optimization only achieves 73.9% sensitivity, 75.0% specificity and 74.4% accuracy. We select the optimal feature subset of the data through MARL and optimize SVM by GA. More remarkably, this method is compared with other dimension reduction ones, such as PCA and the encoder of the VAE. Specifically, the SVM optimized by GA and MARL reached the 100% specificity, 91.3% sensitivity, 94.9% accuracy and 0.870 cross-validation score, achieving the best performance. In summary, the proposed method can improve the NPSLE diagnosis efficiently. In the future, we plan to apply this method to practical medical diagnosis. In the application scenario, we will evaluate the risk of information systems using this technology to prevent some incidents and ensure patient safety ^[49]. Since this technology involves critical infrastructure (healthcare), it is necessary to be regulated.

However, our study still has some limits. On the one hand, further studies and comparisons are needed in patients with HCs, non-NPSLE and NPSLE. Previous studies have found that the cerebral perfusion patterns of NPSLE and non NPSLE patients are different ^[50]. Besides, NPSLE is significantly associated with CII and depression, fatigue and pain ^[51], while non-NPSLE patients and controls don't have it. On the other hand, we do not have sufficient data. We know the importance of abundant labeled data for obtaining good fitting of classifier models ^[52], but the data collection process is complex, and the patients are limited. To address this issue, we use stratified sampling to minimize sampling error and improve the model's generalization ability. By contrast, K-fold cross-validation does not consider the problem of label distribution. Furthermore, to effectively use limited data and overcome overfitting, we use the stratified sampling strategy five times so that the evaluation results can be as close as possible to the model's performance on the test set. The results show that the proposed model can perform well on the unknown data.

In conclusion, our model can perform well in the condition of limited samples with missing values of metabolic features, proving its robustness. Therefore, the proposed algorithm provides a new perspective for early diagnosis and intervention for the NPSLE and reduces the deviation caused by manual screening with more space for its potential development. For the next step, we prepare to collect more samples and related data, including clinical case data and its corresponding medical images, to develop a multimodal machine learning model.

Acknowledgments

This work was supported by the NSFC (No. 61902232), the 2022 Guangdong Basic and Applied Basic Research Foundation grant to Dr. Teng Zhou, the STU Incubation Project for the Research of Digital Humanities and New Liberal Arts (No. 2021DH-3), the STU Scientific Research Foundation for Talents (No. NTF18006), the 2020 Li Ka Shing Foundation Cross-Disciplinary Research Grant (No. 2020LKSFG05D), and the Education Science Planning Project of Guangdong Province (No. 2018GXJK048).

Conflict of interest

The authors declare that they have no conflict of interest.