Extreme Learning Machine Based Prediction of Soil Shear Strength: A Sensitivity Analysis Using Monte Carlo Simulations and Feature Backward Elimination

: Machine Learning (ML) has been applied widely in solving a lot of real-world problems. However, this approach is very sensitive to the selection of input variables for modeling and simulation. In this study, the main objective is to analyze the sensitivity of an advanced ML method, namely the Extreme Learning Machine (ELM) algorithm under di ﬀ erent feature selection scenarios for prediction of shear strength of soil. Feature backward elimination supported by Monte Carlo simulations was applied to evaluate the importance of factors used for the modeling. A database constructed from 538 samples collected from Long Phu 1 power plant project was used for analysis. Well-known statistical indicators, such as the correlation coe ﬃ cient (R), root mean squared error (RMSE), and mean absolute error (MAE), were utilized to evaluate the performance of the ELM algorithm. In each elimination step, the majority vote based on six elimination indicators was selected to decide the variable to be excluded. A number of 30,000 simulations were conducted to ﬁnd out the most relevant variables in predicting the shear strength of soil using ELM. The results show that the performance of ELM is good but very di ﬀ erent under di ﬀ erent combinations of input factors. The moisture content, liquid limit, and plastic limit were found as the most critical variables for the prediction of shear strength of soil using the ML model.


Introduction
In the design phase of various large-scale construction projects (highways, roads, high rise buildings) and geotechnical structures (earth dams, retaining walls), shear strength is an important factor used to define the capability of soil foundations [1]. The shear strength of soil is determined using Mohr-Coulomb criteria through two parameters, namely unit cohesion (c) and internal friction angle (ϕ) in the case of normal soil or only unit cohesion (c) in the case of sandy soil [1]. However, to determine these parameters, the consuming time and costly experiments are often carried out in namely Extreme Learning Machine (ELM) algorithm for prediction of soil shear strength. The main contribution of this study to the knowledge body is that (i) it proposes a soft computing technique (ELM) for quick and accurate prediction of soil shear strength considering more input parameters, which is limited or not easy to be done by using the empirical correlation equation, (ii) it evaluates for the first time the performance of ELM under different combination of input parameters using Monte Carlo simulation and Feature Backward Elimination, which will help in suitable selection of parameters for prediction of soil shear strength using soft computing techniques. For this aim, data of 538 soil samples collected from the Long Phu 1 power plant project, Long Phu district, Soc Trang province, Vietnam were used for generating the datasets used in the modeling. Well-known statistical indicators, such as the correlation coefficient (R), root mean squared error (RMSE), and mean absolute error (MAE), were utilized to evaluate the performance of the ELM algorithm under sensitivity analysis.

Methodology
In order to address this problem, the methodology of the present study contains several main steps such as (1) construction of the database: Input parameters, namely the clay content, moisture content, specific gravity, void ratio, liquid limit, and plastic limit were gathered from technical reports. The considered output variable of this work is the shear strength of soil, (2) ELM algorithm was firstly optimized by an analysis concerning the number of neurons used in the model, (3) after the optimal number of neurons of ELM successfully found, it was used to perform the backward elimination in combination with Monte Carlo simulation, (4) using six types of criteria, namely the maximum value of R, minimum values of RMSE and MAE, average values of R, RMSE and MAE, the elimination of input variables was decided by majority vote.

Data Collection and Preparation
Data used in this study were collected from the Long Phu 1 power plant project (longitude of 9 • 59 07.3"N and latitude of 106 • 04 48"E) located at the southern side of the Hau river, Long Duc commune, Long Phu district, Soc Trang province, Vietnam ( Figure 1). Union of Engineering Geology, Construction and Environment (UGCE) was in charge of the soil investigation works. In addition, a program of the additional soil investigation was carried out by UGCE in April and August 2011, including exploratory borings, field testing, and soil laboratory testing to provide the information relating to the soil conditions of foundation design and construction of the project, and these data were extracted to generate the datasets for the modeling of soil shear strength prediction in this study.
Datasets of 538 soil samples were extracted from the project and used in this study. In datasets, variables, such as moisture content (%), clay content (%), void ratio, plastic limit (%), liquid limit (%), and specific gravity, were used as inputs, and shear strength was used as output. Table 1 shows the initial statistical analysis of the dataset, including the unit and coding of each variable. It is seen that statistical information, such as average, standard deviation, and quantiles of all variables, is fully exposed. For illustration purposes, Figure 2 presents the corresponding histograms of all variables used in this study, as well as the scatter plot between input variables and output response. It can be observed that the distribution of clay covered a wide range between 0 and 65 (mm), the liquid limit from 20 to 65 (%), and the plastic limit from 15 to 35 (%) with a high concentration around 20 (%). Most of the specific gravity values were in the 2.6-2.7 range, whereas the void ratio covered between the 0.5-1.0 range and a low concentration of values was around the 1.75 range. It can also be observed that there is no direct relationship between inputs and output response. Thus, it can be stated that the choice of variables in this study is relevant and suitable [12].  Table 1 shows the initial statistical analysis of the dataset, including the unit and coding of each variable. It is seen that statistical information, such as average, standard deviation, and quantiles of all variables, is fully exposed. For illustration purposes, Figure 2 presents the corresponding histograms of all variables used in this study, as well as the scatter plot between input variables and output response. It can be observed that the distribution of clay covered a wide range between 0 and 65 (mm), the liquid limit from 20 to 65 (%), and the plastic limit from 15 to 35 (%) with a high concentration around 20 (%). Most of the specific gravity values were in the 2.6-2.7 range, whereas the void ratio covered between the 0.5-1.0 range and a low concentration of values was around the 1.75 range. It can also be observed that there is no direct relationship between inputs and output response. Thus, it can be stated that the choice of variables in this study is relevant and suitable [12].
In order to validate the efficiency of the developed ML model, a sub-dataset calling testing part was made, exhibiting 30% (161 samples) of the total 538 configurations. It is worth noticing that such a rate of testing/training was recommended in the literature when developing ML-based models [13][14][15][16][17]. On the other hand, in order to reduce fluctuations within the dataset in training the ML model, as the variables have different ranges of values, all variables were scaled into the range of [0, 1] in order to avoid an unexpected jump in optimizing weight parameters of the models [13,[18][19][20]. The scaling process of a variable x is expressed by Equation (1), and it involves two parameters, α and β, as indicated in Table 1. Precisely, α is the minimum value of the dataset and β is the maximum value.
original scaled In order to validate the efficiency of the developed ML model, a sub-dataset calling testing part was made, exhibiting 30% (161 samples) of the total 538 configurations. It is worth noticing that such a rate of testing/training was recommended in the literature when developing ML-based models [13][14][15][16][17]. On the other hand, in order to reduce fluctuations within the dataset in training the ML model, as the variables have different ranges of values, all variables were scaled into the range of [0, 1] in order to avoid an unexpected jump in optimizing weight parameters of the models [13,[18][19][20]. The scaling process of a variable x is expressed by Equation (1), and it involves two parameters, α and β, as indicated in Table 1. Precisely, α is the minimum value of the dataset and β is the maximum value.    Histograms of the parameters used in this study and correlation graphs with the output: (a,b) Clay, (c,d) moisture content, (e,f) specific gravity, (g,h) void ratio, (i,j) liquid limit, (k,l) plastic limit, and (m) shear strength of soil.

Extreme ML-Based Modeling
Extreme Learning Machine (ELM) is a single hidden layer feedforward neural network (SLFN). The performance of SLFN should be suitable for the system, which can be modeled for data such as critical value, weight, activation function. Therefore, higher learning can be done. In gradient-based learning approaches, all of these parameters are reiteratively modified for each appropriate value. Unlike feedforward neural networks (FNN), which are renewed based on the gradient, in the ELM process, the output weights are analytically built while the input weights are randomly chosen. For an analytic learning process, success rate increases thanks to a strong reduction of the resolution time and the error value. ELM can be introduced to choose a linear function for activating cells in a hidden layer, maybe use non-linear (such as sigmoid and sinusoidal), non-derivatized, or intermittent activation functions [21,22]. ELM algorithm can be shown in the following equations: ( ) kG/cm 2 Figure 2. Histograms of the parameters used in this study and correlation graphs with the output: (a,b) Clay, (c,d) moisture content, (e,f) specific gravity, (g,h) void ratio, (i,j) liquid limit, (k,l) plastic limit, and (m) shear strength of soil.

Extreme ML-Based Modeling
Extreme Learning Machine (ELM) is a single hidden layer feedforward neural network (SLFN). The performance of SLFN should be suitable for the system, which can be modeled for data such as critical value, weight, activation function. Therefore, higher learning can be done. In gradient-based learning approaches, all of these parameters are reiteratively modified for each appropriate value. Unlike feedforward neural networks (FNN), which are renewed based on the gradient, in the ELM process, the output weights are analytically built while the input weights are randomly chosen. For an analytic learning process, success rate increases thanks to a strong reduction of the resolution time and the error value. ELM can be introduced to choose a linear function for activating cells in a hidden layer, maybe use non-linear (such as sigmoid and sinusoidal), non-derivatized, or intermittent activation functions [21,22]. ELM algorithm can be shown in the following equations: g(w n,1 x n + a 1 ) · · · g(w n,m x m + a m ) where α i is the weights between the input layer and the hidden layer and α j is the weights between the output layer and the hidden layer, a j is the critical value of the neurons in the hidden layer, g(.) activation function. Input layer weights (w i,j ) and bias (a j ) are randomly selected. At the beginning of the input layer neuron number (n) and hidden-layer neuron number (m), the activation function (g(.)) is selected. To construct the ELM algorithm, the database was split into a training dataset (70% data) and the remaining data (30%) for building and validation of the ELM model.

Backward Elimination-Based Sensitivity Analysis
Backward elimination, belonging to the wrapper methods, is basically the opposite of the forward selection approach [23]. Precisely, all input variables are firstly chosen, then the most unimportant of the variables are removed one by one in this case [24]. For strategic choices of the process, relative importance of an input variable can be obtained by eliminating an input variable and assessing the influence on the model to be retrained without it or by examining the effect of each input variable on the output by the sensitivity analysis method. In the filtering strategies, the least relevant candidates will be deleted repeatedly until the optimal criteria are satisfied. The process of backward elimination can be summarized in Figure 3.
where αi is the weights between the input layer and the hidden layer and αj is the weights between the output layer and the hidden layer, aj is the critical value of the neurons in the hidden layer, g(.) activation function. Input layer weights (wi,j) and bias (aj) are randomly selected. At the beginning of the input layer neuron number (n) and hidden-layer neuron number (m), the activation function (g(.)) is selected. To construct the ELM algorithm, the database was split into a training dataset (70% data) and the remaining data (30%) for building and validation of the ELM model.

Backward Elimination-Based Sensitivity Analysis
Backward elimination, belonging to the wrapper methods, is basically the opposite of the forward selection approach [23]. Precisely, all input variables are firstly chosen, then the most unimportant of the variables are removed one by one in this case [24]. For strategic choices of the process, relative importance of an input variable can be obtained by eliminating an input variable and assessing the influence on the model to be retrained without it or by examining the effect of each input variable on the output by the sensitivity analysis method. In the filtering strategies, the least relevant candidates will be deleted repeatedly until the optimal criteria are satisfied. The process of backward elimination can be summarized in Figure 3.

Monte Carlo Simulations
Monte Carlo method is one of the most widely used techniques for propagating the input variability on the output results [25][26][27][28][29]. Regarding, for instance, the field of geotechnical engineering, Pham et al. [11] applied the Monte Carlo method for accounting variability of various content properties of soil on the prediction of its mechanical behavior under compression during a highway project. In another attempt for steel structures, the Monte Carlo technique was employed by Le et al. [15] in order to quantify the robustness of hybrid ML models for predicting the critical buckling load of structural members. For typical construction and building materials, such as concrete, many studies involving Monte Carlo technique were introduced in the literature, taking into account the variability in the input space. For instance, Wang et al. [30] quantified the size effect of random aggregates and pores on the mechanical properties of concrete. Jaskulski et al. [31] proposed a probabilistic analysis for concrete subjected to shear. So far, numerical prediction models involving Monte Carlo method could strongly explain the variation of the output results through statistical analysis.
Monte Carlo method is extremely robust and efficient for calculating the propagation of the input variability on the output results, especially using ML models [11,32]. The main idea of the Monte Carlo method is to repeat realizations randomly in the input space and then calculate the corresponding output through the simulation model [33,34]. Therefore, this numerical technique exhibits a high ability in parallel computing [35][36][37][38]. A concept of using the Monte Carlo method is presented in Figure 4, involving a two-dimensional input space with a typical probability distribution. Monte Carlo method is one of the most widely used techniques for propagating the input variability on the output results [25][26][27][28][29]. Regarding, for instance, the field of geotechnical engineering, Pham et al. [11] applied the Monte Carlo method for accounting variability of various content properties of soil on the prediction of its mechanical behavior under compression during a highway project. In another attempt for steel structures, the Monte Carlo technique was employed by Le et al. [15] in order to quantify the robustness of hybrid ML models for predicting the critical buckling load of structural members. For typical construction and building materials, such as concrete, many studies involving Monte Carlo technique were introduced in the literature, taking into account the variability in the input space. For instance, Wang et al. [30] quantified the size effect of random aggregates and pores on the mechanical properties of concrete. Jaskulski et al. [31] proposed a probabilistic analysis for concrete subjected to shear. So far, numerical prediction models involving Monte Carlo method could strongly explain the variation of the output results through statistical analysis.
Monte Carlo method is extremely robust and efficient for calculating the propagation of the input variability on the output results, especially using ML models [11,32]. The main idea of the Monte Carlo method is to repeat realizations randomly in the input space and then calculate the corresponding output through the simulation model [33,34]. Therefore, this numerical technique exhibits a high ability in parallel computing [35][36][37][38]. A concept of using the Monte Carlo method is presented in Figure 4, involving a two-dimensional input space with a typical probability distribution. In this work, the statistical convergence of Monte Carlo simulations has been investigated using the following equation [18,32,39,40]: In this work, the statistical convergence of Monte Carlo simulations has been investigated using the following equation [18,32,39,40]: where G is the mean value of the considered random variable G and n MC is the number of Monte Carlo runs. This convergence function provides efficient information related to the computational time, reliability results for further statistical analysis.

Performance Evaluation
To validate the predictive capability of the models, Mean Absolute Error (MAE), the Pearson correlation coefficient (R), and Root Mean Squared Error (RMSE) were selected and used, as these validation criteria are popular in evaluating the ML models. Basically, R indicates the statistical relationship between the actual values of experiments and the predicted values of the models [41]. Its absolute values range from 0 to 1 where 0 shows an inaccurate correct model and 1 indicates an accurate model. Higher R values indicate better performance of the models. RMSE indicates the average squared difference between the actual and predicted values [42]. In the case of MAE, it shows the average of absolute difference between predicted and actual values [43]. In general, RMSE and MAE show the error evaluation of the models. Thus, lower RMSE and MAE values indicate better performance of the models. Calculation of these values (R, RMSE, and MAE) can be carried out using the following equations: where m is defined as the number of samples, r i and t are defined as the values and means of the predicted shear strength, respectively, and t i and t are the values and mean of the actual shear strength, respectively.

Validation of ELM with Various Number of Neurons
Validation of ELM was conducted by performing 1000 simulations to each of 12 ELM architectures, where the number of neurons varied from 5 to 60 with a step of five neurons. Overall, the total number of simulations was 12,000, taking into account the random sampling index of the dataset. The results with respect to R, RMSE, and MAE are plotted in Figure 5, where the red squares represent the average values, and the blue bars show the standard deviation with respect to 1000 simulations. On the basis of average values and standard deviation of R, RMSE, and MAE, it is found that the optimal number of neurons is in the range of 15 to 25. The best performance of ELM is with 20 neurons, where the highest value of R and lowest values of RMSE and MAE were obtained. Moreover, the standard deviation of ELM using 20 neurons over 1000 random simulations is also smaller compared to other neuron options.
The obtained values of average and standard deviation of RMSE are 0.1082 and 0.0231, whereas those of MAE are 0.0857, 0.0231, respectively, and those of R are 0.9218 and 0.0167, respectively. Overall, the performance of ELM is good for the prediction of shear strength of soil, and the number of 20 neurons used for training ELM was selected as an optimal choice for further investigations.

Sensitivity Analysis Using Backward Elimination and Monte Carlo Simulations
The sensitivity analysis by performing backward elimination with the help of Monte Carlo simulations is carried out in this section. A number of four scenarios (Scenarios 1 to 4), corresponding to each input space after the elimination process, was defined. The "Scenario 0" refers to the case using the initial input space without excluding any variables. The "Scenario 1" consisted of six different input spaces containing only five variables, in which each variable was excluded from the corresponding input space. For instance, six input spaces considered in this case were: (i) X2, X3, X4, X5, X6; (ii) X1, X3, X4, X5, X6; (iii) X1, X2, X4, X5, X6; (iv) X1, X2, X3, X5, X6; (v) X1, X2, X3, X4, X6; (vi) X1, X2, X3, X4, X5. Similarly, the "Scenario 2", "Scenario 3", and "Scenario 4" corresponded to the cases with five, four, three input spaces, respectively (Figure 3). The summarized input space and the four scenarios could be illustrated in Figure 3. The following sections are dedicated to each step of the backward elimination process.
3.2.1. Reduction of the Input Space from 6 to 5 Variables (Scenario 1) The first step of backward elimination consists of quantifying the performance of ELM in predicting the shear strength of soil by excluding each variable successively in the input space of the database. Thus, a number of 6000 simulations (six input spaces x 1000 simulations) were performed in excluding successively from input X1 to X6. The results are plotted in Figure 6 for average values of R, RMSE, and MAE (red squares), standard deviation (blue bars) and min, max values (orange bars). Detailed values with respect to six elimination indicators are summarized in Table 2. For the sake of comparison, the discontinuous black lines represent the corresponding values of the criteria for the case of using all input variables (Scenario 0). On the basis of average values of R, it is observed

Sensitivity Analysis Using Backward Elimination and Monte Carlo Simulations
The sensitivity analysis by performing backward elimination with the help of Monte Carlo simulations is carried out in this section. A number of four scenarios (Scenarios 1 to 4), corresponding to each input space after the elimination process, was defined. The "Scenario 0" refers to the case using the initial input space without excluding any variables. The "Scenario 1" consisted of six different input spaces containing only five variables, in which each variable was excluded from the corresponding input space. For instance, six input spaces considered in this case were: (i) X 2 , X 3 , X 4 , X 5 , X 6 ; (ii) X 1 , X 3 , X 4 , X 5 , X 6 ; (iii) X 1 , X 2 , X 4 , X 5 , X 6 ; (iv) X 1 , X 2 , X 3 , X 5 , X 6 ; (v) X 1 , X 2 , X 3 , X 4 , X 6 ; (vi) X 1 , X 2 , X 3 , X 4 , X 5 . Similarly, the "Scenario 2", "Scenario 3", and "Scenario 4" corresponded to the cases with five, four, three input spaces, respectively (Figure 3). The summarized input space and the four scenarios could be illustrated in Figure 3. The following sections are dedicated to each step of the backward elimination process.
3.2.1. Reduction of the Input Space from 6 to 5 Variables (Scenario 1) The first step of backward elimination consists of quantifying the performance of ELM in predicting the shear strength of soil by excluding each variable successively in the input space of the database. Thus, a number of 6000 simulations (six input spaces x 1000 simulations) were performed in excluding successively from input X 1 to X 6 . The results are plotted in Figure 6 for average values of R, RMSE, and MAE (red squares), standard deviation (blue bars) and min, max values (orange bars). Detailed values with respect to six elimination indicators are summarized in Table 2. For the sake of comparison, the discontinuous black lines represent the corresponding values of the criteria for the case of using all input variables (Scenario 0). On the basis of average values of R, it is observed that the performance of the ELM algorithm in excluding clay content (X 1 ) slightly decreased from 0.9218 (simulation with six inputs) to 0.9203 (simulation with five inputs except for clay content). For the remaining cases (excluding from X 2 to X 6 ), the performance of ELM decreased more significantly. Similar observations were noticed taking the average values of RMSE and MAE. Indeed, it was found that excluding clay content (X 1 ) reduced the ELM prediction performance with RMSE decrease from 0.1082 to 0.0925, and MAE decreased from 0.0857 to 0.0722 while comparing the cases of all input variables and without clay content in the input space. Besides, taking the maximum values of R or minimum value of MAE as an indicator, plastic limit (X 6 ) was the variable to be excluded. However, taking the minimum value of RMSE as an indicator, the specific gravity (X 3 ) was the variable to be excluded. Finally, the elimination decision was made based on the majority vote between indicators, where clay content (X 1 ) was selected to be a less important variable compared with other variables for predicting soil shear strength.     The second step of backward elimination consists of the assessment of ELM capability in predicting the shear strength of soil by excluding each of the remaining inputs (X2 to X6). Thus, a number of 5000 simulations (5 input spaces × 1000 simulations) were performed. The results of average and standard deviation values of R, RMSE and MAE are displayed in Figure 7. Detailed values with respect to six indicators are summarized in Table 3. The discontinuous black lines represent the error criteria values for the case without using clay content (X1) as input variable. On the basis of average values of R, it is observed that the performance of the ELM algorithm in excluding void ratio (X4) decreased from 0.9203 (simulation with five inputs except clay content) to 0.9188 (simulation with four inputs except clay content and void ratio). For the remaining cases (excluding X2, X3, X5, and X6), the performance of ELM exhibited lower values (Table 3). Similar remarks were observed for the average values of RMSE and MAE. Indeed, it was found that excluding the void ratio made inconsiderable changes with RMSE (increase from 0.0925 to 0.0957) and MAE (increase from 0.0722 to 0.0751) with respect to the cases of all input variables without clay content as a variable. Interestingly, taking the maximum values of R, or minimum values of RMSE and MAE as indicators, the void ratio was also the variable to be excluded. The elimination at this stage revealed that void ratio (X4) is a less important variable compared with other variables (X2, X3, X5 and X6) in predicting the soil shear strength.  The second step of backward elimination consists of the assessment of ELM capability in predicting the shear strength of soil by excluding each of the remaining inputs (X 2 to X 6 ). Thus, a number of 5000 simulations (5 input spaces × 1000 simulations) were performed. The results of average and standard deviation values of R, RMSE and MAE are displayed in Figure 7. Detailed values with respect to six indicators are summarized in Table 3. The discontinuous black lines represent the error criteria values for the case without using clay content (X 1 ) as input variable. On the basis of average values of R, it is observed that the performance of the ELM algorithm in excluding void ratio (X 4 ) decreased from 0.9203 (simulation with five inputs except clay content) to 0.9188 (simulation with four inputs except clay content and void ratio). For the remaining cases (excluding X 2 , X 3 , X 5 , and X 6 ), the performance of ELM exhibited lower values (Table 3). Similar remarks were observed for the average values of RMSE and MAE. Indeed, it was found that excluding the void ratio made inconsiderable changes with RMSE (increase from 0.0925 to 0.0957) and MAE (increase from 0.0722 to 0.0751) with respect to the cases of all input variables without clay content as a variable. Interestingly, taking the maximum values of R, or minimum values of RMSE and MAE as indicators, the void ratio was also the variable to be excluded. The elimination at this stage revealed that void ratio (X 4 ) is a less important variable compared with other variables (X 2 , X 3 , X 5 and X 6 ) in predicting the soil shear strength.    The third step of backward elimination consists of predicting the shear strength of soil by successively excluding the remaining inputs (X 2 , X 3 , X 5 , and X 6 ). This induces a total number of 4000 simulations (4 input spaces × 1000 simulations) to be performed. The results of average and standard deviation values of R, RMSE, and MAE are plotted in Figure 8. Detailed simulation results with respect to six indicators are summarized in Table 4. The discontinuous black lines represent the error criteria values for the simulation without using clay content (X 1 ) and void ratio (X 4 ) as input variables. With respect to the average values of R, it is observed that the performance of ELM algorithm in excluding plastic limit (X 6 ) slightly decreased from 0.9188 (simulation with four inputs without clay content and void ratio) to 0.9164 (simulation with three inputs except for clay content, void ratio, and plastic limit). On the contrary, different remarks were observed for the average values of RMSE and MAE. Precisely, it was found that excluding specific gravity made inconsiderable changes with RMSE (increase from 0.0931 to 0.0993) and MAE (increase from 0.0732 to 0.0778). More importantly, taking the maximum values of R, or minimum values of RMSE and MAE as indicators, specific gravity was the variable need to be eliminated. The backward elimination at this stage revealed that the specific gravity (X 3 ) is less important input variable compared with other variables (X 2 , X 5 , and X 6 ) in predicting the shear strength of soil. The third step of backward elimination consists of predicting the shear strength of soil by successively excluding the remaining inputs (X2, X3, X5, and X6). This induces a total number of 4000 simulations (4 input spaces × 1000 simulations) to be performed. The results of average and standard deviation values of R, RMSE, and MAE are plotted in Figure 8. Detailed simulation results with respect to six indicators are summarized in Table 4. The discontinuous black lines represent the error criteria values for the simulation without using clay content (X1) and void ratio (X4) as input variables. With respect to the average values of R, it is observed that the performance of ELM algorithm in excluding plastic limit (X6) slightly decreased from 0.9188 (simulation with four inputs without clay content and void ratio) to 0.9164 (simulation with three inputs except for clay content, void ratio, and plastic limit). On the contrary, different remarks were observed for the average values of RMSE and MAE. Precisely, it was found that excluding specific gravity made inconsiderable changes with RMSE (increase from 0.0931 to 0.0993) and MAE (increase from 0.0732 to 0.0778). More importantly, taking the maximum values of R, or minimum values of RMSE and MAE as indicators, specific gravity was the variable need to be eliminated. The backward elimination at this stage revealed that the specific gravity (X3) is less important input variable compared with other variables (X2, X5, and X6) in predicting the shear strength of soil.    The final step of backward elimination in this study consists of performing the prediction by excluding one of the remaining inputs (X2, X5, and X6). At this stage, the total number of 3000 simulations (3 input spaces × 1000 simulations) was performed. The results of average and standard deviation values of R, RMSE, and MAE are plotted in Figure 9. Detailed simulation results with respect to six indicators are summarized in Table 5. The discontinuous black lines represent the error criteria values for the simulation without using clay content (X1), specific gravity (X3) and void ratio (X4) as input variables. With respect to the average values of R, it is observed that the performance of ELM algorithm in excluding plastic limit (X6) slightly decreased from 0.9138 (simulation with four inputs without clay content, specific gravity, and void ratio) to 0.8815 (simulation with two inputs moisture content and liquid limit). On the contrary, different remarks were observed for the average  The final step of backward elimination in this study consists of performing the prediction by excluding one of the remaining inputs (X 2 , X 5 , and X 6 ). At this stage, the total number of 3000 simulations (3 input spaces × 1000 simulations) was performed. The results of average and standard deviation values of R, RMSE, and MAE are plotted in Figure 9. Detailed simulation results with respect to six indicators are summarized in Table 5. The discontinuous black lines represent the error criteria values for the simulation without using clay content (X 1 ), specific gravity (X 3 ) and void ratio (X 4 ) as input variables. With respect to the average values of R, it is observed that the performance of ELM algorithm in excluding plastic limit (X 6 ) slightly decreased from 0.9138 (simulation with four inputs without clay content, specific gravity, and void ratio) to 0.8815 (simulation with two inputs moisture content and liquid limit). On the contrary, different remarks were observed for the average values of RMSE and MAE. Precisely, it was found that excluding the liquid limit made inconsiderable changes with RMSE (increase from 0.0993 to 0.1334) and MAE (increase from 0.0778 to 0.0778). On the other hand, taking the maximum values of R, or minimum values of RMSE and MAE as indicators, the plastic limit was the variable need to be eliminated. Overall, it could be considered that the plastic limit is less important than the liquid limit and moisture content in predicting the shear strength of soil, and thus, it can be concluded that moisture content is the most important factor for the prediction of the shear strength of soil.

Performance of ELM in Predicting the Shear Strength of Soil
Overall, the performance of the ELM algorithm in predicting the shear strength of soil is satisfactory. From Scenarios 1 to 4, in reducing the input space from six to three variables, the accuracy of ELM was reasonably accepted. Precisely, the average values of R decreased from 0.9218 to 0.8815, those of RMSE were varied in the range of 0.1082-0.1334 and an increase of MAE values from 0.0857 to 0.1093 between Scenario 0 and Scenario 4, respectively. ELM algorithm could thus be considered as a good predictor to deal with the soil shear strength problem. Moreover, the computation time of ELM is very fast in comparison with other common ML methods such as ANN, ANFIS, or SVM [11]. For illustration purposes, one simulation using ELM in this study took only less than 0.1 seconds, which could be very efficient for massively parallel computing. The total 30,000 simulations in this study were conducted in just a few hours in an Intel Xeon E3-1505M V5 2.80GHz computer using eight processors. The reason lies in the main concept of ELM, by random initialization of the single hidden layer feedforward NN for the weights and biases [44].

Reliability of the Predicted Results by Monte Carlo Approach
The reliability of an ML algorithm, represented by the convergence of simulation results, is crucial for any analysis. Presented in the previous sections, the convergence analysis via Monte Carlo simulation is performed in order to provide additional information on the prediction capability of the ELM model. It is worth noticing that, for the sake of simplicity, only results of R and RMSE are shown in Figure 10. Overall, the R values over five simulation scenarios were converged after only 200 runs (in a 1% range compared with the average values of R), whereas smaller fluctuation can be achieved over 700 simulations in all the cases. For RMSE values, the fluctuations were observed in the 4% range compared with the corresponding average values of RMSE. More stable results were achieved when the number of simulations exceeded 800. Interestingly, in excluding the liquid limit (X 5 ), the plastic limit (X 6 ) or moisture content (X 2 ) from the input space seemed to highly increase the fluctuations of the statistical analysis. This could be another confirmation to strengthen the conclusion of the backward selection feature in this study, as these three variables were considered as important. The fluctuation analysis of results was in very good agreement with the convergence analysis. Thus, performing backward elimination coupling with Monte Carlo simulation as a support decision indicator, an in-depth point of view on the importance of variables could be revealed.

Backward Elimination Criteria-Based Sensitivity Analysis
The backward elimination process, belonging to the wrapper methods, has been empirically proven that it obtains subsets with better performance than certain feature selection methods (i.e., filters methods) as the obtained subsets are evaluated by real modeling algorithms [24]. Moreover, comparing with forward selection, backward elimination finds a stronger subset of features because of the full assessment of features during the selection process [23]. Last but not least, backward elimination could reveal information not only on the importance but also the unimportance of input variables. However, selecting the criteria to remove such variables from the input space is crucial. In this study, six criteria were proposed and independently evaluated to provide useful information to the final decision. The average values of R, RMSE, or MAE could be potential candidates as they reflect the global, converged performance of the ELM algorithm over a sufficient number of simulations. It is found that, in all cases, the average values of RMSE and MAE were the reliable indicators for the elimination process. The average value of R was failed in one case (Scenario 3), similar to the minimum values of RMSE and MAE (Scenario 1).
The results show that the values of R attained negative values in several cases. This could be the reason for several sudden changes in the statistical convergence curves (Figure 9), making the average value of R not a good indicator. This is quite familiar for all ML algorithms performed with Monte Carlo simulation, as the values of R could be stable (ANFIS model in [45], SVM algorithm in [11]) or unstable (ANN algorithm in [45]). The values of RMSE and MAE were noticed to be more stable, in which extreme or outlier values were not observed. This could be because using the average values

Backward Elimination Criteria-Based Sensitivity Analysis
The backward elimination process, belonging to the wrapper methods, has been empirically proven that it obtains subsets with better performance than certain feature selection methods (i.e., filters methods) as the obtained subsets are evaluated by real modeling algorithms [24]. Moreover, comparing with forward selection, backward elimination finds a stronger subset of features because of the full assessment of features during the selection process [23]. Last but not least, backward elimination could reveal information not only on the importance but also the unimportance of input variables. However, selecting the criteria to remove such variables from the input space is crucial. In this study, six criteria were proposed and independently evaluated to provide useful information to the final decision. The average values of R, RMSE, or MAE could be potential candidates as they reflect the global, converged performance of the ELM algorithm over a sufficient number of simulations. It is found that, in all cases, the average values of RMSE and MAE were the reliable indicators for the elimination process. The average value of R was failed in one case (Scenario 3), similar to the minimum values of RMSE and MAE (Scenario 1).
The results show that the values of R attained negative values in several cases. This could be the reason for several sudden changes in the statistical convergence curves (Figure 9), making the average value of R not a good indicator. This is quite familiar for all ML algorithms performed with Monte Carlo simulation, as the values of R could be stable (ANFIS model in [45], SVM algorithm in [11]) or unstable (ANN algorithm in [45]). The values of RMSE and MAE were noticed to be more stable, in which extreme or outlier values were not observed. This could be because using the average values of RMSE, MAE is better than that of R. Interestingly, the maximum values of R were also found as a reliable elimination indicator in this study. In fact, choosing ML model with the highest accuracy (maximum value of R or R 2 ) to perform further investigation is very common in the literature, for instance, in [17]. The maximum values of R or minimum values of RMSE and MAE only represented the best performance of the ELM algorithm over a certain number of Monte Carlo simulations. Even though the number of 1000 runs in each case was proven to be statistically satisfied, these values could change when increasing the number of simulations. Therefore, the use of maximum values of R or minimum values of RMSE and MAE as backward elimination indicators still needs further investigations.

Importance of Input Factors for Prediction of Soil Shear Strength
Validation of the importance of input factors for developing and applying ML models is an essential task which will help in selecting the most suitable factors used for more effective and accurate modeling and prediction. In this study, with the help of a combination approach of EML, Monte Carlo and backward elimination, the importance of input factors for prediction of soil shear strength was validated and determined (Table 6). Except for Scenario 2, the study found that the moisture content, liquid limit and plastic limit were the three most important input variables. Considering Scenarios 1 and 2, the liquid limit is more important than others, however, Scenarios 3 and 4 demonstrated that the moisture content had a stronger effect on predicting the soil shear strength than others. In general, water-related factors are the most important parameters affecting the soil shear strength and the performance of the predictive ML model. It is reasonable because water can significantly reduce the friction and link between soil particles, thus, the shear strength of less-water soil will be higher than those of more-water soil. The finding of this study is also comparative with other published studies [46,47]. Table 6. Order of importance over four scenarios in this study. Importance  1  2  3  4  5  6 Scenario 1 X 1 X 4 X 3 X 6 X 2 X 5 Scenario 2 X 4 X 2 X 3 X 6 X 5 -Scenario 3 X 3 X 6 X 5 X 2 --Scenario 4 X 5 X 6 X 2 ---

Order of
The limitation of this study is that the feature interaction might be a phenomenon that slightly changed the order of importance in the present study [48]. It occurs when such feature has a little correlation with the predicted target but highly correlated when treating with another feature, so that excluding these types of features could reduce the performance of the ML algorithms [48]. With respect to the change of ML algorithms due to the elimination of one input, the readers could refer to the literature [14].

Conclusions
In this study, the sensitivity of an advanced ML method, namely the ELM algorithm under different feature selection scenarios for prediction of shear strength of soil was carried out. Feature backward elimination and Monte Carlo simulations were applied to evaluate the importance of factors used for the modeling. A database including input variables (moisture content (%), clay content (%), void ratio, plastic limit (%), liquid limit (%), and specific gravity) and output variable (shear strength of soil) constructed from 538 samples collected from Long Phu 1 power plant project was used for analysis. Well-known statistical indicators such as R, RMSE, and MAE were utilized to evaluate the performance of ELM algorithm. In each elimination step, the majority vote was selected to decide the variable to be excluded.
The results show that the performance of ELM is good but very different under different combinations of input factors for the prediction of shear strength of soil. The moisture content, liquid limit, and plastic limit were found as the most important variables, and other factors are less important for prediction of shear strength of soil using the ML model. This study might help to select the suitable factors for more quickly and accurately prediction of shear strength of soil using ML models.