Unifying Variable Importance Scores from Different Machine Learning Models Using Simulated Annealing

2.1 Variable importance analysis

It is widely known that a predictive model resulting from a machine learning algorithm is more difficult to interpret than a classical statistical model [7]. The complication of the model comes from the fact that the model does not have an explicit mathematical function mapping a set of predictor variables to a single value of the response variable. The model tends to have a non-linear relationship among variables.

Due to that circumstance, a follow-up analysis is needed to reveal the information from the model [8]. One of the analyses useful to reveal more information from the machine learning models is the variable importance analysis (VIA) [7]. The VIA tries to reveal each predictor variable's contribution in determining the model's performance to predict the output of new observations [9]. The predictor variable with a higher contribution would have a higher ranking among the predictors [8].

Several VIA methodologies can be found in the literature. In general, the VIA procedures produce a set of scores from which the analyst can identify the relative importance of each predictor. Some existing methods are permutation variable importance, Shapley Value, and SHAP. The permutation variable importance [10] has been a methodology that is widely implemented in many situations. This approach could be applied in models generated by any machine learning algorithm such as NN, RF, XGB, and SVM. The other advantage of this methodology is that it requires a reasonably fast computational time.

Suppose we have a predictive model F, which utilizes p predictor variables X₁, X₂, …, X_p to predict a variable response Y. Further, suppose that the matrix of the predictor variable values is X and the vector of variable response is y. The performance of the model F could be measured by a loss function L, which is obtained by comparing the actual response vector y and the prediction response vector $\widehat{\boldsymbol{y}}={\boldsymbol{F(X)}}$. The variable importance score of a predictor variable X_j is obtained by the following steps [9]:

1. Calculate $L^0=\mathcal{L}(\widehat{\boldsymbol{y}}, \boldsymbol{X}, \boldsymbol{y})$; this is the value of the loss function in the original data.
2. Create a matrix $\boldsymbol{X}^{* j}$ by permutating the j^th column of X.
3. Calculate the predictive model of $\widehat{\boldsymbol{y}}^{* j}$ based on $\boldsymbol{X}^{* j}$, which is $\widehat{\boldsymbol{y}}^{* j}=F\left(\boldsymbol{X}^{* j}\right)$.
4. Calculate the second loss function value by comparing y and $\widehat{\boldsymbol{y}}^{* j}: L^{* j}=\mathcal{L}\left(\widehat{\boldsymbol{y}}^{* j}, \boldsymbol{X}, \boldsymbol{y}\right)$.
5. Calculate the importance score of $\boldsymbol{X}^j$ by $s_j=L^{* j}-L^0$ or $s_j=\frac{L^{* j}}{L^0}$.

In a classification problem, the loss function might be 1 – AUC or misclassification rate. A variable with a high value of s_j indicates that X_j is more important because the exclusion of this variable from the model will increase the loss function so that the model performance is low.

Unfortunately, the result of s_j score may vary when we change the machine learning algorithm to produce model F. Even if the dataset used is identical, the importance rankings might differ between those obtained from the random forest methodology and those obtained from other machine learning algorithms. This circumstance could raise an issue in deciding which subset of variables should be considered significant. This paper will discuss the proposed approach to ease analysts' decision-making.

2.2 Notations

Suppose we ran k different machine learning algorithms to generate k predictive models. For each model, we implemented the permutation variable importance procedure to have k sets of variable importance rankings R_i, i = 1, …, k. The set R_i consists of p values or R_i = {s_i1, s_i2, …, s_ip} where s_ij is the rank of the variable importance for predictor variable X_j, j = 1, …, p. Note that s_ij should be a whole number between 1 to p or s_ij Î {1, 2, …, p} and s_ij ≠ s_ik for all j ≠ k.

Recall that instead of the variable importance score, the value of s_ij is the rank of the importance score descendingly. A variable with a smaller s_ij is considered the more important predictor since it is associated with a higher value of variable importance score.

Our proposed method aims to obtain a single set of variable importance ranks from those k sets. In other words, we would like to unify the k sets into one. Let us use the notation of W = {w₁, …, w_p} as the result of the unifying process. As w_j is the importance rank of predictor variable X_j, it also has properties like what s_ij has. The value of w_j should satisfy the conditions

w_j Î {1, 2, …, p} and w_j ≠ w_k for all j ≠ k.

The basic idea of the unification is finding w_j; j = 1, …, p so that the correlation between W and R_i is as high as possible for all i = 1, …, k. In other words, we would like to find a new variable importance rank W to agree highly with all original sets of importance rank R_i’s. Suppose that z_i is the correlation between W and R_i, z_i = cor(W, R_i) for i = 1, …, k. The proposed method tries to reach a situation so that z_i is high for all i = 1, .., k. To ensure that, we define

Z = min {z_i}        (1)

and use it as the objective function of the maximization problem.

The complete formulation of the optimization problem is then could be written as follows:

Find W = {w₁, …, w_p} that maximize Z

Subject to constraints:

w_j Î {1, 2, …, p}; for all j = 1, …, p

w_j ≠ w_k for all j ≠ k

The above optimization problem could be seen as a combinatorial problem since the solution is the permutation of p whole number from 1 to p. Therefore, we propose implementing a simulated annealing algorithm as a metaheuristic approach to handle this optimization problem. The following subsection describes the details of the algorithm.

2.3 Computational procedure

As mentioned, a simulated annealing algorithm is utilized to find the solution for the optimization problem to unify the importance ranks. The basic algorithm is described as follows.

At the initiation step, we define the values of T_max and T_min. Both values determine the number of iterations to improve the solution. Then, an initial solution is generated. In our method, a random permutation of whole numbers 1 to p will be improved during the process. Initially, we define the value T = T_max [11].

The iteration starts by generating a new solution slightly different from w. It can be generated by exchanging the values of two entries in w and named as w’. The positions of the changed entries were randomly chosen. The process continues by replacing w by w’ whenever the performance of w’ is better than w [12, 13]. In this paper, the performance of the solution is Z, the minimum correlation between w and R_i. The higher the Z value, the better the solution.

The reason behind the choice of maximizing the minimum value of correlations is as follows. A good solution of feature ranking is one that has a high correlation with every single result from machine learning algorithms. If the minimum correlation value is large, then all correlation values are larger. Therefore, maximizing minimum correlation means that the algorithm would yield high values for all correlation coefficients.

The second condition to replace w by w’ is a random process. Even though Z(w’) is lower than Z(w), there is a possibility that w’ replaces w with the probability is inversely proportional to the number of iterations and the difference between the performance. It means that if the iteration has gone quite long, the probability of replacing the current solution with the worse one is lower. Also, if the performance of the alternative solution is much worse, then the probability is much lower.

As the iteration goes, the value of T is decaying. The algorithm will stop if the value of T reaches T_min [14], and the most recent w is the final solution containing the set of the importance ranks.

The pseudocode of the algorithm can be written as follows [15]:

Initiation:

set T_min and T_max

set $\theta$, set k

w=w₀ /* Generate initial solution */

Evaluate the solution f(w)

T=T_max /* Initial temperature */      

Repeat

Generate a new solution w’

$\Delta E=Z\left(w^{\prime}\right)-Z(w)$

If $\Delta E>0$ Then s=s’/* Accept new solution */

Else Accept s’ with probability $p=exp \left(-\frac{\Delta E}{k T}\right)$

T=$\theta$T /* New temperature */

Until T < T_min

Output: The best solution is fulfilled

s = current solution; s₀ = initial solution; k = Boltzman constant; T = temperature; s’ = new solution; f(s) = objective value for solution s.

4.1 Empirical data

The data comes from the 2020 national economic survey (Susenas) in West Java, Indonesia. The number of predictor variables is 24, and one response variable has two classes: 0 indicates a food-secure family, and 1 indicates a food-insecure family. The number of observations is 24,679 families. The names of the predictors can be seen in Table 3.

4.2 Variable importance analysis for empirical data

Permutations variable importance (PVI) methods are applied to FIES data using four models and obtained variable importance measures in Figure 2. PVI is prioritized because of simplicity and ease of interpretation. The PVI is conducted through a hyperparameter tunning technique to obtain optimal results [21]. RF, XGB, NN, and SVM are applied because the models are suitable for classification data, the models can handle relationship non-linear between predictor variables and response variable, and they are used widely because of their effectiveness. The order of the variable importance from RF machine learning is X1, which is the most important; X2, the sixth most important; and so on. For XGB variable importance (VI), X1 is in first place, X2 is in second place, and so on, which differs from the results of RF VI. For NN, the ranking of the importance of variables is the same as those of RF except for variables X6, X8, X11, X13, X14, X16, X18, and X19. Most of the VI order from SVM differs from the order of VI from other VI. There are differences in variable importance measure (VIM) in each machine learning model, as in the case of the SHAP-FI analysis in previous studies [5]. This VIM distinction makes it difficult to interpret. Therefore, a merger method is proposed in this study.

Table 3. Predictors, symbol, and measure in Food Insecure Experience Scale (FIES) data

The model's performance is shown in terms of accuracy [22]. The values of accuracy and Area Under the ROC Curve (AUC) of each ML model are quite high. The accuracy values for RF, XGB, NN, and SVM are 0.7465, 0.7826, 0.7465, and 0.7488, respectively, while the AUC is 0.813, 0.858, 0.813, and 0.812. The accuracy values exceed 0.74, while the AUC values exceed 0.81 (Figure 3).

4.3 Implementation of the proposed method uses empirical data

Objective values have converged after the 300th iteration. The objective values are minimum (correlation [s_i', RF-VI]), correlation [s_i', XGB-VI], correlation [s_i', NN-VI], correlation [s_i', SVM-VI]). Correlation [s’, RF-VI] means the correlation coefficient between the proposed variable importance on iteration i and the variable importance of RF. Python software is used to carry out the simulation annealing technique. If the objective value converges, the optimal combined variable importance has been obtained (Figure 4).

2.png

Figure 2. Score variable importance RF, XGB, NN, and SVM

3.png

Figure 3. Accuracy and AUC score of RF, XGB, NN, SVM

The order of the variable importance for each machine learning method produces the joint VIM. The joint variable importance showed a strong Spearman correlation with the machine learning models RF, XGB, NN, and SVM, scoring 0.953; 0.923; 0.936; and 0.926, respectively. The five most important predictor variables in influencing food insecure families are X3 (education of the family head), X1 (house size), X4 (No of my family members have a savings account), X5 (Floor-type of the house), and X7 (Ownership of land). (Table 4).

4.png

Figure 4. Iteration to obtain the proposed optimal combined variable importance measure

Table 4. Machine-learning and simulated annealing (SA) variable importance

4.4 Evaluation of the proposed method with empirical data

The joint variable importance scores can be seen to describe the order of influence of variables. The joint variable importance score is calculated by averaging the scores of the four PVI (Figure 5). Based on the average, X1, X2, X3, X4, X5, X6, X7, X8, and X9 have scores greater than 0.03. A score of 0.03 means the distance between the loss function score and the original data loss function score is 0.03. In other words, the predictors significantly influence the response variable. It works like the elbow method in cluster analysis [23]. If predictors. Next, characteristics of VI-SA are identified through the number of predictors.

A boxplot assesses the SA algorithm's stability in creating VI (Figure 6). The object value (minimum correlation) increased and became more stable when more variables were included. This study looks at a variety of variables, ranging in number from few to numerous. The variables included 5, 10, 15, 20, and 24. In conclusion, various factors can improve the accuracy of the joint variable's important measure. This study examines several different variables.

5.png

Figure 5. Scores of average variable importance of random forest, XGBoost, Neural Network, Support Vector Machine

6.png

Figure 6. Comparison of objective value for the number of variables m = 5, m = 10, m = 15, m = 20, and m = 24 with 50 replications

Evaluation of uncertainty was carried out by observing the SA process 100 times. The predictors have a rating uncertainty. The results of the SA process can produce insignificant different VI. Patterns of median ranks, ranks in percentile 5, and ranks in percentile 95. Ranks of the proposed variable importance have a level of uncertainty. The ranks can different between line p0.05 and p0.95 (Figure 7)

The proposed variable importance measure is seen through individual classification tree algorithms. The classification tree is a method to predict output [24]. The accuracy of the joint VIM increases when the number of variables increases. The accuracy is quite good when the number of predictors exceeds 10 (Figure 8). The accuracy shows that the proposed method is suitable.

7.png

Figure 7. Patterns for ranks of each predictor with 100 replications

8.png

Figure 8. Accuracy scores of simulated annealing variable importance performances use a classification tree

In summary, results of its application to empirical data support that it can be applied, and results of evaluation is satisfied. The machine learning and joined variable importance is optimal when predictors are independent each others. The constituent variable importances should have high accuracies so the proposed variable importance has high accuracy. The optimal accuracy is gotten when the number of predictors exceeds 10 (ten). The proposed variable importance have uncertainty characteristics.