Hybrid mRMR and multi-objective particle swarm feature selection methods and application to metabolomics of traditional Chinese medicine

mRMR

mRMR (Peng, Long & Ding, 2005) is to find the set of features in the original set of features that have the highest relevance to the final output result (Max-Relevance) and the set of features with the minor correlation between features (Min-Redundancy). Assuming that x and y are two random variables, the mutual information is defined as:

(1) $$I(x;y)=\int \int p(x,y)\log {\displaystyle \frac{p(x,y)}{p(x)p(y)}}dxdy.$$

The maximum class correlation $\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{D}\left(\mathrm{S},\mathrm{C}\right)\right)$ for the target class C and a subset of features with S features is defined by the average of the correlations between the selected features $x_i\left(i=1,2,\dots ,S\right)$ and the target class C:

(2) $$max(D(S,C))D={\displaystyle \frac{1}{s}\sum _{x_i\in S}I(x_i,C).}$$

The following equation can evaluate the redundancy of S:

(3) $$minR(s)={\displaystyle \frac{1}{\left|S^2\right|}\sum _{(x_i,x_j)\in S^2}I(x_i,x_j).}$$

The mRMR ranks features by simultaneously maximizing relevance and minimizing redundancy and is expressed in the form of the following equation:

(4) $$\mathrm{m}\mathrm{R}\mathrm{M}\mathrm{R}=max\left[D={\displaystyle \frac{1}{s}\sum _{x_i\in S}I(x_i,C)-R={\displaystyle \frac{1}{\left|S^2\right|}\sum _{(x_i,x_j)\in S^2}I(x_i,x_j)}}\right].$$

Multi-objective optimization

Various complicated tasks need to optimize two or more objectives simultaneously, and the objective functions conflict with each other and cannot explicitly balance them; that is, the optimization of one objective may cause the decay of another objective and cannot make each objective function optimal. The challenge of simultaneously optimizing conflicting objectives within a defined range is called a multi-objective optimization problem (MOOP) (Lücken Von, Brizuela & Barán, 2019). The goal is to find solutions to achieve the best possible trade-off among the competing objectives. When solving practical problems with multi-objectives, feature selection will be considered a maximization or minimization problem. If the number of objectives is n, a standard multi-objective optimization mathematical model should be denoted as shown in Eq. (1). Multi-objective optimization methods aim to obtain a set of non-dominated optimal solutions or a subset of characteristics.

(5) $$\begin{array}{rl}& min(max)F(x)=\{f_1(x),f_2(x)...f_n(x)\}\\ & s.t.x\in \mathrm{\Omega }\end{array}.$$

The Pareto dominance method (Newman, 2005) is a standard method for solving MOOP. It is defined by examining two decision vectors $a,b\in X$, $a$ Pareto dominance $b$, noted as $a>b$, when and only when the following equation, which can also be called $a$ dominates $b$. If no other decision variable can dominate it, then the decision variable is termed a non-dominated solution.

(6) $$\{\mathrm{\forall }i\in \{1,2,...,n\}{f}_i(a)\le f_i(b)\}\wedge \left\{\mathrm{\exists }j\in \left\{1,2,...,n\right\}{f}_j(a)<f_j(b)\right\}.$$

Particle swarm optimization

Particle swarm optimization (PSO), as described in the work previously published by Figueiredo, Ludermir & Bastos-Filho (2016), is a type of swarm intelligence algorithm inspired by the foraging behavior of bird flocks. A weightless particle simulates birds within a flock, possessing two properties: velocity (V) and position (X). The velocity attribute of the particle corresponds to the magnitude or speed at which it moves, while the position attribute indicates the direction or orientation of its movement. Each particle conducts an individual search within the search space to find its personal best solution, which is known as the $P_{best}$. This $P_{best}$ is then shared among all particles in the swarm. Together, they determine the current global best ( $G_{best}$) solution, representing the overall best solution found by the entire swarm. Each particle adjusts its velocity and position based on its $P_{best}$ and the current $G_{best}$ solution shared among the particles. The position of the particle is $X_i=(x_{i1},x_{i2},...,x_{iD})$, the velocity is $v_i=(v_{i1},v_{i2},...,v_{iD})$, the individual extremum is $P_i=(p_{i1},p_{i2},...,p_{iD})$ and the global optimal solution is $P_g=(p_{g1},p_{g2},...,p_{gD})$. Equations give the updated velocity and position formulas:

(7) $$V_{id}(t+1)=w\ast V_{id}(t)+c_1\ast r_1\ast (p_{id}(t)-x_{id}(t))+c_2\ast r_2\ast (p_{gd}(t)-x_{id}(t))$$

(8) $$x_{id}(t+1)=x_{id}(t)+V_{id}(t+1)i=1,2,...,N,d=1,2,...,D$$

where t denotes the number of iterations and N is the size of the population.

Fitness function

Currently, the PSO algorithm is mainly applied to continuous optimization problems, and the MOOP has been studied less. Feature selection has two conflicting objectives, feature subset minimization, and performance maximization (Zhang et al., 2023), and is a typical MOOP.

In this study, the number of feature subsets and the root mean square error (RMSE) of the evaluation index of the regression task are chosen as the objective functions of the MOPSO, where the number of feature subsets is denoted as $f_1$. The RMSE is denoted as $f_2$, as shown in Eqs. (9) and (10).

(9) $$f_1=num$$

(10) $$f_2=\sqrt{{\displaystyle \frac{1}{n}{\sum }_{i=1}^n{\left(Y-\hat{Y}\right)}^2}}$$where, $num$ denotes the number of selected features, $Y$ is the predicted value, and $\hat{Y}$ is the actual value.

Multi-objective particle swarm optimization based on dynamic acceleration factor and nonlinear decreasing inertia weights (CMOPSO)

To address the problem of the MOPSO algorithm quickly falling into the local optimum, this article improves the convergence of MOPSO by using the dynamic acceleration factor and nonlinear decreasing inertia weights to make the locally optimal particles jump out. The specific flow of the algorithm is outlined in Algorithm 1.

The equations for the nonlinear decreasing inertia factor (w) value and the dynamic acceleration factor are shown in Eqs. (11) and (12), where r1 and r2 refer to two random values within [0, 1], t denotes the number of iterations, and T denotes the total number of iterations. The change in the magnitude of the nonlinearly decreasing inertia factor value (w) with an increasing number of iterations is shown in Fig. 2. In the literature (Feng, Chen & Guo, 2006), the range of the acceleration factor is obtained: when $c_1=2.75\sim 1.25$ and $c_2=0.5\sim 2.25$, the effect is better, so the parameters $c_{1f}=2.75$, $c_{1i}=1.25$, $c_{2f}=0.5$ and $c_{2i}=2.25$ are set. Then, the size changes with the number of iterations, as shown in Fig. 3, with the number of iterations increasing linearly and decreasing, with the number of iterations increasing and growing linearly.

(11) $$w=(w_{max}-w_{min})\ast (t/T-1{)}^2+w_{min}$$

(12) $$\begin{array}{rl}& c_1=c_{1f}+(c_{1i}-c_{1f})\ast (t/T)\\ & c_2=c_{2f}+(c_{2f}-c_{2i})\ast (t/T)\end{array}$$

Hybrid mRMR and multi-objective particle swarm feature selection methods

To address the feature selection problem in HDSS settings, this article introduces a two-stage method called MCMOPSO. During the first phase of MCMOPSO, the mRMR algorithm measures the correlation between features and the target feature and the redundancy among features. By prioritizing maximum relevance and minimum redundancy, the features are sorted. The Wrapper strategy is then employed to eliminate irrelevant and partially redundant features, enabling the adaptive discovery of the best-performing candidate feature subset. In Phase 2, the remaining redundant features are removed from the MOPSO based on dynamic acceleration factors and nonlinear decreasing inertia weights, and a subset of features with a small number and high accuracy is selected.

Assuming that there are $m$ features, $n$ samples, a subset of features after data pre-processing $F=(f_1,f_2,...,f_m)$, and sample features after pre-processing $Y=(y_1,y_2,...,y_n{)}^T$, the specific construction process of the MCMOPSO model is as follows:

Phase 1: Filter coupled with Wrapper to eliminate irrelevant and partially redundant features

Step1. mRMR calculates correlation and redundancy: the mRMR scores between each feature $f_i$ and the target features $Y$ and $f_i$ are calculated to obtain the score sequence $mRM{R}_{list}=(f_1:mRMR(f_1,Y),f_2:mRMR(f_2,Y),...,f_m:mRMR(f_m,Y))$, and the score sequence $mRM{R}_{list}$ is sorted in descending order.

Step2. Determine the subset of candidate features: To address the challenge of determining the number of features and threshold values in the Filter method; this article incorporates the Wrapper approach to determine the optimal number of features to retain adaptively. In this method, a certain number of features are added by a forward search strategy. Then, PLS is used to model the regression of feature subsets, and finally the one with the smallest RMSE value is the candidate feature subset $F_{can}$.

Phase 2: Use a heuristic search algorithm to remove redundancy again

Step3. Eliminate redundant features: The proposed MOPSO optimization algorithm, which incorporates dynamic acceleration factors and non-linearly decreasing inertia weights, is employed to remove redundant features from the candidate feature subset ( $F_{can}$). This process results in the identification of the optimal feature subset ( $F_{best}$).

Step4. Model evaluation: The effectiveness and timeliness of the model are evaluated by assessing the performance of the optimal feature subset.

The complete MCMOPSO algorithm is shown in Algorithm 2:

Validation of CMOPSO’s validity on regular datasets

CMOPSO changes the acceleration factor and inertia weights basis on MOPSO, using dynamic acceleration factor and nonlinear decreasing inertia weights to make better diversity for particles to jump out of the local optimum in the early stage and speed up the particles in the later stage to improve the performance of MOPSO and make it have a better optimization finding accuracy. The regular dataset on UCI in Table 2 was used to verify the validity of CMOPSO and compared with two feature selection methods, approximate Markov blanket (AMB) and DA2MBL1, in the literature (Li et al., 2022). The DA2MBL1 method obtains similar feature groups by approximate Markov blanket clustering of features, imposes L1 regular term constraints on each similar feature group, and combines the coordinate descent method of solving to make the regression coefficients of the redundant features compressed to 0 to achieve the purpose of deleting redundant features. To avoid the chance of CMOPSO, each data is run ten times, and group’s solution with the smallest $f_1$ or $f_2$ is taken each time. The average of ten times is finally taken, and the result after 10-fold cross-validation is shown in Table 7. Bolded text is the optimal value of the result.

Table 7, shows that the CMOSPO algorithm has better results on the SPmath, SPPort, and Rbuild datasets. Especially on the Rbuild dataset, y1 and y2 are left with only three features from 103 features after the CMOPSO algorithm removes irrelevant and redundant features, and the RMSE of y1 decreases from the original 589.1403 to 207.4673, and the RMSE of y2 decreases from the original 55.0247 to 34.8772. From the three datasets in Table 7, it can be seen that the CMOPSO algorithm has the lowest number of features and the lowest RMSE among the three algorithms, so CMOPSO may have better results in removing irrelevant and redundant features on regular datasets compared to AMB and DA2MBL1 algorithms.

Parameter setting

This article will involve two multi-objective optimization algorithms: MOPSO and Non-dominated Sorting Genetic Algorithm-II (NSGA-II) (Deb et al., 2002). NSGA-II is one of the popular multi-objective methods that uses a fast, non-dominated sorting algorithm and introduces an elitist strategy. To demonstrate the effectiveness of the MCMOPSO proposed in this article, the unmodified MOPSO and NSGA-II were used for experimental comparison. The number of iterations set for fairness is 300, the particle swarm size is 100, and the remaining parameters are set as shown in Table 8. The conventional parameters of MOPSO and CMOPSO are the same.

Validation of MCMOPSO effectiveness on high dimensional small sample datasets

To address the feature selection problem for metabolomics high-dimensional small-sample data, this article proposes the MCMOPSO method using a two-stage feature selection framework, tested on two metabolomics data and 10 conventional high-dimensional small-sample data. In this section, two multi-objective methods, NSGA-II and MOPSO, are used to compare with the second-stage MCMOPSO. To better compare the three algorithms, all use mRMR to adaptively select the candidate feature set as input features and the number of feature subsets and RMSE as two optimization objectives. To avoid the chance of the algorithm results, each comparison algorithm is run independently ten times on each dataset. The best Pareto solution with the same number of features in the results of the ten runs is averaged, and the smaller the RMSE is, the better it is for the same number of features. The final results are presented in a bar chart, and the error line is added to indicate the difference between the results of the ten runs and the average. This time, the error metric in the graph is standard deviation, and the shorter length of the error line indicates that the results of each run are very close to the mean and the algorithm has good stability.

Figure 4 shows four endogenous substance data to verify the validity of MCMOPSO. The horizontal coordinate of this bar graph indicates the number of feature subsets, and the vertical coordinate indicates the RMSE size. As the three algorithms perform better when the number of features in the Endo-y1 data is taken as 9 and 12, the average RMSE of NSGA-II is 331.53 and 380.02, the average RMSE of MCMOPSO is 302.51 and 314.65, respectively, and the average RMSE of MOPSO is 306.85 and 330.95, respectively. For NSGA-II, MCMOPSO and MOPSO have error values of 23.085, 11.78, and 0 for 9 features and 0, 0, and 10.34 for 12 features, respectively. To compare of the mean RMSE, the MCMOPSO algorithm performs the best, and the NSGA-II algorithm performs the worst for both features. As for the error line, each algorithm’s performance varies for different features. CMOPSO algorithm has a low error line for nine features, and the MOPSO algorithm has a low error line for 12 features. However, the NSGA-II algorithm has a higher error line for both features.

As can be observed from the results in Fig. 4, the average RMSE of MCMOPSO on the four endogenous substances data is smaller. Still, when Endo-y1 retains eight ions, the RMSE is higher than that of NSGA-II but slightly lower than that of MOPSO with a similar magnitude of error. In the Endo-y1 data, MCMOPSO's average RMSE is better than the two comparison algorithms, but the error is more significant when 11 ions are selected. MCMOPSO in the Endo-y3 data outperforms the comparison algorithms in the average RMSE and the error value. The average RMSE is lower than NSGA-II on the Endo-y4 data except for selecting 13 ions, which is lower than the average RMSE of NSGA-II and MOPSO.

In addition, in this article, we conducted independent t-test for the MCMOPSO algorithm against the NSGA-II algorithm and the MOPSO algorithm, respectively, and obtained the corresponding p-values, denoted as p1 and p2, respectively, which are shown in the upper right corner of each figure in Fig. 4. In addition to performing hypothesis tests, we obtained two other evaluation metrics through regression analysis: the R-square and the MAE. We averaged the solutions obtained from one run and ten runs for each dataset and computed the final average of their R-squared and MAE. The results for R-squared and MAE are tabulated for comparison in File S1 (Table S5). From the p-value results for datasets y1, y2, and y3, our algorithm shows a significant advantage (p-value < 0.05) under certain comparison conditions. While the R-square for all four datasets is larger than the two comparison algorithms, the MAE for y1, y3, and y4 is smaller than the two.

In summary, the MCMOPSO algorithm outperforms the two comparison algorithms on endogenous substances overall, but the significance difference is less pronounced and only significant on some of the datasets.

Figure 5 shows the results of the runs on the four pharmacodynamic indicators in exogenous substances; from the figure MCMOPSO algorithm on Exo-y1, Exo-y3, and Exo-y4, it can be seen that the average RMSE is lower than the two comparative algorithms and better performance. Still, on Exo-y2, the performance is not as good as the NSGA-II algorithm but better than the MOPSO algorithm. From the error line, it can be seen that MCMOPSO has fewer errors on Exo-y2-y4 but more on Exo-y1, especially when 10 ions are selected. The MCMOPSO algorithm outperforms NSGA-II and MOPSO algorithms on exogenous substance data.

From the p-values obtained during the hypothesis testing, all p-values are less than 0.1 except for p2 on y1, which is greater than 0.1. and p1 on the four datasets of exogenous substances are less than 0.05. From the results of R-squared and MAE, the R-squared of y2, y3, and y4 are more significant than the two comparative algorithms, and the MAE of y2 and y3 are less than the two comparative algorithms. There is a substantial difference between MCMOPSO and NSGA-II in the exogenous substance dataset, but the critical difference with MOPSO on this dataset is not so noticeable. From the results of multiple regression evaluation metrics, it can be seen that MCMOPSO outperforms the two comparative algorithms in most of the datasets and, therefore, outperforms them in overall performance.

Figure 6 shows the experimental results on ten regular HDSS datasets. The MCMOPSO algorithm significantly outperforms the NSGA-II and MOPSO algorithms regarding average RMSE and error on data BlogTe1-3, BlogTe6, and BlogTe10. The BloTe5 algorithm does not perform as well as the two comparative algorithms when taking three ions but outperforms them when taking 4–6 ions. Similarly, on BlogTe9 data the MCMOPSO algorithm has a lower average RMSE than NSGA-II and MOPSO algorithms when the number of ions taken is three, four, and 11. On BlogTe4, BlogTe7, and BlogTe8, the MCMOPSO algorithm does not perform as well as the two comparison algorithms, but the errors are similar.

The p1 value is less than 0.05 for six data on the public dataset, i.e., more than half of the datasets have a significant difference between MCMOPSO and NSGA-II, and the p2 value is less than 0.05 for four datasets. From the results of the R-square and the MAE, it can be seen that the R-square of the MCMOPSO algorithm outperforms the two comparative algorithms for five datasets, and the MAE outperforms the two comparative algorithms for six datasets (Table S6 in File S1). Therefore, the MCMOPSO algorithm is slightly better than the NSGA-II algorithm and MOPSO algorithm in overall performance.

In summary, the MCMOPSO algorithm proposed in this article is slightly better than the NSGA-II and MOPSO algorithms in dealing with HDSS data, and the error is more minor in most cases. However, there is a significant error in some data or when choosing a specific number of features, so the algorithm’s stability needs to be improved. The experimental comparison shows that the difference between NSGA-II, MCMOPSO, and MOPSO is more evident in the Exo dataset, especially the BlogTe dataset. This is because the Endo dataset has higher dimensionality than the Exo and BlogTe datasets, and NSGA-II maintains the diversity of the population by using non-dominated ordering and crowding degree distance. More dimensions provide more possibilities and diversity in higher dimensional space due to larger solution space. So, NSGAII works slightly better with higher dimensional endo datasets. Different algorithms may show different advantages and disadvantages for a MOOP like feature selection. This leads to more apparent differences between NSGA-II, CMOPSP, and MOPSO on Exo datasets, especially BlogTe datasets.

Comparison of algorithm performance

Time efficiency is also an essential element of feature selection research, and this article will compare and analyze the algorithm’s running time and time complexity. Table 9 lists the running times of the proposed algorithm and the two comparison algorithms. The table shows that CMOPSO has much less running time than NSGA-II on regular high-dimensional data and higher-dimensional data in metabolomics, with 13 of them having less running time than MOPSO. Because CMOPSO changes the acceleration factor based on the original MOPSO, it is a little faster and takes less time to run. In terms of time complexity, assuming that the size of the particles in the particle swarm and the size of the population in NSGA-II are both n, and the particle size and initial population size are m, the time complexity of MOPSO and CMOPSO is $max[O(mn),O(n),O(n^2)]$, and the time complexity of NSGA-II is $max[O(mn),O(n^2),O(n^3)]$. Comparing the time complexity, the time complexity of MOPSO is lower than that of NSGA-II. The experimental results show that CMOPSO has better calculation efficiency in the feature selection algorithm based on multi-objective optimization.

To present a more intuitive comparison of the running time of the three algorithms, the time comparison in Table 9 is visualized in the form of a bar chart, and the results are shown in Fig. 7. The horizontal coordinates represent different datasets, the vertical coordinates represent the running time, and the three colors represent the three algorithms. From the figure, it can be seen that the running time of the NSGA-II algorithm is greater than that of CMOPSO and MOPSO algorithms, whereas, on the datasets Endo-y1-y4, Exo-y2-y3, BlogTe1, and BlogTe5-10, the running time of CMOPSO algorithm is less than that of MOPSO algorithm. Therefore, the CMOPSO algorithm is slightly better than the two compared algorithms in terms of time performance.