Stability Risk Assessment of Underground Rock Pillars Using Logistic Model Trees

Pillars are important structural elements that provide temporary or permanent support in underground spaces. Unstable pillars can result in rock sloughing leading to roof collapse, and they can also cause rock burst. Hence, the prediction of underground pillar stability is important. This paper presents a novel application of Logistic Model Trees (LMT) to predict underground pillar stability. Seven parameters—pillar width, pillar height, ratio of pillar width to height, uniaxial compressive strength of rock, average pillar stress, underground depth, and Bord width—are employed to construct LMTs for rock and coal pillars. The LogitBoost algorithm is applied to train on two data sets of rock and coal pillar case histories. The two models are validated with (i) 10-fold cross-validation and with (ii) another set of new case histories. Results suggest that the accuracy of the proposed LMT is the highest among other common machine learning methods previously employed in the literature. Moreover, a sensitivity analysis indicates that the average stress, p, and the ratio of pillar width to height, r, are the most influential parameters for the proposed models.


Introduction
Pillars are important structural elements in rock and coal underground mining since they can provide temporary or permanent support for tunneling and mining work [1][2][3]. To maximize the extraction rate of mineral resources, dimensions (width and height) of underground pillars must be confined to a certain range, which should still assure stability and safe working conditions throughout the entire life of such projects [4]. Unstable pillars can result in rock cracking, causing roof collapse, and they could even trigger rock bursts when large amounts of accumulated elastic energy are released [3,5]. Increases of mining depth can lead to more frequent pillar instability and failure events; it is therefore expected that new design challenges must be addressed so that new methods to assess rock pillar stability be discussed and developed.
The most common parameter in pillar design and stability estimation is the factor of safety (FoS), which can be calculated as the ratio of pillar strength to applied pillar stress. Two main methods -the tributary area theory, and numerical simulations-are often employed for pillar stress calculation [3,6]; whereas several empirical methods have been developed for pillar strength estimation, as listed in Table 1. Three basic 'effects' related with the ratio of pillar height to width-the linear shape effect, size effect, and power shape effect-were summarized by Lunder [6]; similarly, Lunder [6] worked out an empirical equation to assess pillar strength, considering confined and unconfined strength and different ratios of pillar height to width. Similarly, Ahmad and Al-Shayea et al. [7] proposed new coal pillar strength formulae that consider the interface friction between the roof/floor and pillar, which also have a linear and power shape effect. Notation: Ps, pillar strength (MPa); UCS, unconfined compressive strength of a cubic pillar specimen (MPa); w, pillar width (m); h, pillar height (m); σ 1 , major principal stress (MPa); σ 3 , Minor principal stress (MPa); σ c , unconfined compressive strength of intact rock (MPa); m and s, empirically derived constants based on rock mass quality of pillar material; κ, pillar confinement/friction factor; C pav , average σ 3 /σ 1 ratio across the mid-height centerline of the pillar.
Theoretically, a rock or coal pillar is predicted as 'stable' when the value of FoS is greater than 1; otherwise, it is considered 'unstable'. However, such crisp boundaries are often unreliable, as unstable pillars often occurs when the FoS value is above 1 (Lunder [6]), and uncertainties in the 'capacity' (strength or resisting force) and 'demand' (stress or disturbing force) of a pillar lead to the FoS being described as a statistical distribution.
Numerical simulation has also been employed to evaluate pillar stability. For instance, York [15] and Mortazavi and Hassani et al. [16] employed FLAC 2D to estimate the stability of deep-level pillars considering pillar deformations under natural loading conditions. Li and Li et al. [17] employed FLAC 3D to obtain the minimum required thickness of crown pillars under sea water pressure; and Martin and Maybee [18] conducted two-dimensional finite element analyses using Hoek-Brown parameters for pillar strength. Similarly, the finite element method and the discrete fracture network (DFN) approach were integrated by Deng and Yue et al. [19] and Elmo and Stead [20] to study the failure of jointed pillars. However, although such numerical methods are powerful tools to model pillar behavior, they have limitations due to the complexity of the models involved and the need for specific calibration.
Many machine learning algorithms have been applied to predict pillar stability. For instance, Tawadrous and Katsabanis [21] employed Artificial Neural Networks (ANN). Recio-Gordo and Jimenez [22] presented a probabilistic model based on the theory of linear classifiers that can be used to make probabilistic predictions of pillar behavior in longwall and retreat room and pillar mining. Zhou and Li et al. [2] employed Fisher Discriminant Analysis (FDA) and Support Vector Machines (SVM), whereas Wattimena and Kramadibrata et al. [23] and Wattimena [4] applied logistic regression algorithms. Ghasemi and Ataei et al. [24] and Ghasemi and Ataei et al. [25] developed logistic regression and fuzzy logic algorithms in room-and-pillar coal mines. Zhou and Li et al. [3] employed six machine learning algorithms to evaluate the stability of a rock pillar, and Ghasemi and Kalhori et al. [26] compared Decision Tree (J48) and SVM algorithms for hard rock pillar stability assessment. Mohanto and Deb [27] developed a plastic damage index for assessing rib pillar stability in underground metal mines using multi-variate regression and artificial neural network techniques. Ahmad and Al-Shayea et al. [7] employed random trees and C4.5 decision trees for underground pillar prediction, and Liang and Luo et al. [28] used new GBDT, XGBoost, and LightGBM algorithms. Dai and Shan et al. [29] studied an intelligent identification method for coal pillar stability in a fully mechanized caving face of thick coal seams. With a combination of finite difference methods, neural networks, and Monte Carlo simulation techniques, Li and Zhou et al. [30] studied underground mine hard rock pillar stability.
However, despite their reliable and accurate results, most algorithms are not easily applicable in practice due to their complex training and modeling procedures, and to their 'black box' features. The model tree algorithm-which jointly uses decision trees and linear regression methods-was proposed by Quinlan [31] to overcome such limitations. This algorithm has been successfully applied to numerous geotechnical problems, such as UCS and Young's Modulus estimation [32] and hydraulic conductivity [33]. However, one of the key obstacles to using tree algorithms for pillar stability evaluation is that outcomes should be discrete values such as 'stable', 'unstable' or 'failed'. Therefore, those discrete targets must be transformed into continuous values, so that they can be trained using a model tree algorithm for regression analysis.
We propose a Logistic Model Trees (LMT) method [34] to predict rock and coal pillar stability. This algorithm has the advantage of dealing with the classification problem by jointly using a tree model and a logistic regression algorithm, making it a rational choice in classification and decision-making. One main application of the LMT in geotechnical engineering has been landslide susceptibility prediction [35][36][37], but it has not yet been applied to predict underground pillar stability.
The flowchart of this research is shown in Figure 1. The rest of this paper is organized as follows. Section 2 discusses two databases employed in this research, and the selection of input parameters. The LMT method is briefly introduced in Section 3. Section 4 discusses the performance and validation of proposed models, and risk and feature importance analyses.
regression and artificial neural network techniques. Ahmad and Al-Shayea et al. [7] employed random trees and C4.5 decision trees for underground pillar prediction, and Liang and Luo et al. [28] used new GBDT, XGBoost, and LightGBM algorithms. Dai and Shan et al. [29] studied an intelligent identification method for coal pillar stability in a fully mechanized caving face of thick coal seams. With a combination of finite difference methods, neural networks, and Monte Carlo simulation techniques, Li and Zhou et al. [30] studied underground mine hard rock pillar stability.
However, despite their reliable and accurate results, most algorithms are not easily applicable in practice due to their complex training and modeling procedures, and to their 'black box' features. The model tree algorithm-which jointly uses decision trees and linear regression methods-was proposed by Quinlan [31] to overcome such limitations. This algorithm has been successfully applied to numerous geotechnical problems, such as UCS and Young's Modulus estimation [32] and hydraulic conductivity [33]. However, one of the key obstacles to using tree algorithms for pillar stability evaluation is that outcomes should be discrete values such as 'stable', 'unstable' or 'failed'. Therefore, those discrete targets must be transformed into continuous values, so that they can be trained using a model tree algorithm for regression analysis.
We propose a Logistic Model Trees (LMT) method [34] to predict rock and coal pillar stability. This algorithm has the advantage of dealing with the classification problem by jointly using a tree model and a logistic regression algorithm, making it a rational choice in classification and decision-making. One main application of the LMT in geotechnical engineering has been landslide susceptibility prediction [35][36][37], but it has not yet been applied to predict underground pillar stability.
The flowchart of this research is shown in Figure 1. The rest of this paper is organized as follows. Section 2 discusses two databases employed in this research, and the selection of input parameters. The LMT method is briefly introduced in Section 3. Section 4 discusses the performance and validation of proposed models, and risk and feature importance analyses.

Database Description
We compiled two databases that include recorded rock and coal pillar stability events from different types of underground mines. Database A was developed using observations collected by Lunder [6]. It includes 178 case histories (60 stable, 50 unstable, and 68 failed cases) that came from six hard rock mines in Canada, South Africa, and Sweden. Five main parameters-pillar width, w; pillar height, h; ratio of pillar width to height, r; uniaxial compressive strength of rock, UCS; and average pillar stress, p-that could affect rock pillar stability were selected and analyzed.
Similarly, Database B was developed based on Van der Merwe [38] and Van der Merwe [39], who provided information on 351 case histories of coal pillars (including 274

Database Description
We compiled two databases that include recorded rock and coal pillar stability events from different types of underground mines. Database A was developed using observations collected by Lunder [6]. It includes 178 case histories (60 stable, 50 unstable, and 68 failed cases) that came from six hard rock mines in Canada, South Africa, and Sweden. Five main parameters-pillar width, w; pillar height, h; ratio of pillar width to height, r; uniaxial compressive strength of rock, UCS; and average pillar stress, p-that could affect rock pillar stability were selected and analyzed.
Similarly, Database B was developed based on Van der Merwe [38] and Van der Merwe [39], who provided information on 351 case histories of coal pillars (including 274 stable pillars and 77 failed cases) in South African coal mines. It originally contains four parameters related to coal pillar stability analyses -underground depth, H; pillar width, w; pillar height (mining height), h; bord width, B-so the ratio of pillar width to height (r) has been computed to be employed as a parameter in this database.
Although there are other features that could affect the rock and coal pillar stability, such as the geology structure conditions, water and gas contents, other stresses, and fracture information from CT or SEM images [14,[40][41][42], that information is not available for the current pillar stability analysis. Since our models are trained using a limited number of data points are empirical methods, it is expected that their predictions are improved when a more extensive data set and/or more features are employed. Thus, a better calibration of model parameters, and more elaborate models could be established.
The boxplots of both databases are shown in Figures 2 and 3, and their statistical features are listed in Tables 2 and 3; solid black spots in Figures 2 and 3 denote "outliers". The first and third quartiles are indicated using horizontal lines at the bottom and top of the boxes, whereas the bold lines inside the boxes represent median values. Similarly, pillars with 'failed', 'stable', and 'unstable' conditions are shown separately, and labelled using F, S, and U, respectively. when a more extensive data set and/or more features are employed. Thus, a better calibration of model parameters, and more elaborate models could be established.
The boxplots of both databases are shown in Figure 2; Figure 3, and their statistical features are listed in Table 2; Table 3; solid black spots in Figure 2; Figure 3 denote "outliers". The first and third quartiles are indicated using horizontal lines at the bottom and top of the boxes, whereas the bold lines inside the boxes represent median values. Similarly, pillars with 'failed', 'stable', and 'unstable' conditions are shown separately, and labelled using F, S, and U, respectively.      Based on the review and analysis of common techniques to estimate pillar stability using FoS, two categories of parameters were employed in our analyses. One category is the "strength" related parameters: pillar width, w; pillar height or mining height, h; ratio of pillar width to height, r; and uniaxial compressive strength of rock, UCS. The other category includes parameters related to the applied "stress," including average pillar stress, p; underground depth, H; and bord width, B. With the availability of parameters for the two databases, Models 1 (w, h, r, UCS and p) and 2 (w, h, H, B and r) were established. The detailed analyses of these parameters are discussed in the following subsections.

Strength-Related Parameters
Pillar width w and pillar height h are two key geometrical parameters for stability analysis. The empirical approaches (linear, power, or shape forms) in Table 1 all included w, h, and r for pillar strength calculation. Zhou and Li et al. [2] employed w, h, r, and UCS as inputs for FDA and SVM models; Wattimena and Kramadibrata et al. [23] and Wattimena [4] used r as one of the parameters for regression analysis. Moreover, pillar shape, Based on the review and analysis of common techniques to estimate pillar stability using FoS, two categories of parameters were employed in our analyses. One category is the "strength" related parameters: pillar width, w; pillar height or mining height, h; ratio of pillar width to height, r; and uniaxial compressive strength of rock, UCS. The other category includes parameters related to the applied "stress," including average pillar stress, p; underground depth, H; and bord width, B. With the availability of parameters for the two databases, Models 1 (w, h, r, UCS and p) and 2 (w, h, H, B and r) were established. The detailed analyses of these parameters are discussed in the following subsections.

Strength-Related Parameters
Pillar width w and pillar height h are two key geometrical parameters for stability analysis. The empirical approaches (linear, power, or shape forms) in Table 1 all included w, h, and r for pillar strength calculation. Zhou and Li et al. [2] employed w, h, r, and UCS as inputs for FDA and SVM models; Wattimena and Kramadibrata et al. [23] and Wattimena [4] used r as one of the parameters for regression analysis. Moreover, pillar shape, as represented by the r parameter, could have an influence on increased strength [3]. Therefore, we consider pillar width w, pillar height h, the ratio r, and UCS as "strength"related parameters for Model 1. As Database B does not include UCS values, only w, h, and r are employed for Model 2.

Stress-Related Parameters
The actual stress acting on a pillar may depend on multiple factors such as in situ stresses, mining induced stress changes, geological features, shape and orientation of pillars, spatial relationship between pillars and mine openings, and groundwater conditions [6]. The tributary area theory introduced by Bunting [43] calculates the average stress acting on a pillar by simply supposing that the pillar supports its "share" of applied load. Salamon and Munro [11] proposed a more simplified equation: where γ is the unit weight of the rock mass, H is underground depth, B is bord width between pillars, and w is pillar width. For Model 1, the average stress has been calculated using the tributary area theory and two-dimensional or three-dimensional finite element modeling in the database, so that it can be directly applied for prediction. For Model 2, H, B, and w are employed to compute the stress condition for coal pillar stability. Table 4 lists the selected parameters for each model, their corresponding data availability, and types. Theoretically, there may be additional indicators or parameters, such as geological condition, pillar orientation, or mining methods, that could have an influence on pillar stability; however, collecting these data is a major challenge and they have therefore not been employed in this work. Table 4. Parameter selection for each model.

Parameters
Model 1 Model 2

Logistic Model Trees
Logistic regression is a simple prediction tool, with advantages such as stability, low variance, and time-efficient training [44], but its prediction results are often biased. Decision trees are another machine learning method for searching a less restricted space of candidate models and capturing nonlinear patterns in a database; they exhibit low bias but high variance and instability, and are hence prone to overfitting. Therefore, the Logistic Model Trees (LMT) approach was proposed by Landwehr and Hall et al. [34]. It builds from the Model Tree approach proposed by Quinlan [31] to deal with regression problems with the joint use of linear regression and decision tree models, and is extended to classification problems. A brief introduction to LMT is presented in this section, while a more detailed description can be found in the seminal work of Landwehr and Hall et al. [34].

Tree Structure
A LMT includes a standard decision tree structure with logistic regression functions generated at the leaves. It has a set of inner or non-terminal nodes N and a set of leaves or terminal nodes T. Unlike a common decision tree or model tree, each leaf t of a LMT model has correlative logistic regression (LR) functions, instead of having classification labels or linear regression functions. For instance, suppose that the input vector X is (X 1 , X 2 ) and the output target is Y. The whole instance space is denoted as S, and such space can be divided into several subspaces S t. A simple input space split into seven subspaces is shown in Figure 4. We can note that: linear regression functions. For instance, suppose that the input vector X is (X1, X2) and the output target is Y. The whole instance space is denoted as S, and such space can be divided into several subspaces St. A simple input space split into seven subspaces is shown in Figure 4. We can note that: For the seven divided subspaces in Figure 4, LR functions are determined in the model. The tree structure is shown in Figure 5. For the seven divided subspaces in Figure 4, LR functions are determined in the model. The tree structure is shown in Figure 5.

Logistic Function
Unlike traditional forms of logistic regression, the LogitBoost algorithm for fitting additive logistic regression models proposed by Friedman and Hastie et al. [45] is employed here for model construction. The prediction probability is presented in Equation where G is the output, J are the class labels, X are the inputs, and Fj (x) are the functions to be trained in the leaves of the tree by the LMT, as:

Logistic Function
Unlike traditional forms of logistic regression, the LogitBoost algorithm for fitting additive logistic regression models proposed by Friedman and Hastie et al. [45] is employed here for model construction. The prediction probability is presented in Equation (3).

of 19
where G is the output, J are the class labels, X are the inputs, and F j (x) are the functions to be trained in the leaves of the tree by the LMT, as: where m is the number of iterations, f mj are the functions of input variables, α are the intercepts and coefficients of the linear function, and s are the variables of the subset S t at the leaf t.

LMT Training
A LMT can be established with the following steps: initial tree growing, tree splitting and stopping, and tree pruning. In this section, the basic idea is introduced; the reader is referred to Landwehr and Hall et al. [34] for detailed information.
The M5P method commonly employed for tree growing can build a standard tree first; then, a logistic regression model can be obtained at every node [32,46]. As this approach merely trains the model using case histories at each node in isolation, without considering the surrounding tree structure, another approach-one that can incrementally refine logistic model fit at high levels-is employed, so that the LogitBoost algorithm can iteratively change F j (x) to increase the fit in a natural way [34]. The function f mj is added to F j by altering one of the coefficients in the function or bringing in another variable (see Equation (4)).
Therefore, first, a LR tree is built in the root using proper iteration numbers determined by cross-validation in the initial growing process. Next, using the C4.5 splitting law [31] to raise the purity of the classification variable, the tree begins to grow by resembling advisable subsets (t) from database (S) to the children nodes. In the children nodes, the logistic regression functions are built by running the LogitBoost algorithm with the consideration of logistic model, weights, and probability estimate performed in the last iteration at the parent node. Then, another splitting process is performed. Considering reliable model fitting, the tree will stop splitting when a node has less than 15 cases. After the tree is built, the tree pruning method is employed to trade off tree size and to reduce model complexity, while still maintaining predictive accuracy. After different pruning scheme experiments, Landwehr and Hall et al. [34] employed the CAERT pruning method [47] to make pruning decisions considering training error and model complexity. Following these three steps, logistic model trees can be established.

Development of LMT Models for Pillar Stability Prediction
We employed the WEKA (Waikato Environment for Knowledge Analysis) [48] software to build models using the two databases collected. Since some case histories in Database A are incomplete (see Table 2), 153 case histories-55 stable cases, 34 unstable cases, and 64 failed cases-were employed as the training set, and nine case histories not used for training were randomly selected for validation. For Database B, 266 stable cases and 73 failed cases were employed as the training set and 12 case histories were selected for the validation work.
The ratio of stable to failed case numbers in the training set of Model 2 is about 3.64, indicating that the distribution of those two classes requires a cost-sensitive approach to overcome this 'lack of balance' problem [49,50]. Therefore, we consider the cost during the training process of Model 2. The cost-sensitive classifier is matched with the LMT algorithm using WEKA. As Model 2 is a binary prediction problem, a 2 × 2 cost matrix, shown in Table 5, is utilized during the training step. The value of C (F|S) is set to 4, so that false positive cases are four times more likely than false negative cases. During tree model training, each terminal node (or leaf) is trained and updated using logistic regression models (see Section 3). Considering predictive performance, easily applicable tree structures, and the total number of training case histories, the minimum number of instances were set to 15 and 40 for Models 1 and 2, respectively. The trees produced using LMT are demonstrated in Figures 6 and 7. Model 1 contains eight logistic functions (LMs), whereas Model 2 includes ten LMs. The detailed expressions are listed in Tables 6 and 7.    As the variance inflation factor (VIF) can measure the correlation and strength of correlation between the variables in a regression model, we used R programing (car package) for the calculation of the VIF value; the results are as follows.
VIF > 10 indicates a potential multicollinearity problem in the dataset [35]. All the VIF values in Table 8 are smaller than 10, which is acceptable for regression analysis.
It should be noted that some functions in Table 6; Table 7 do not include all the parameters selected. For instance, the LM3 function for unstable pillars in Table 6 does not consider UCS. In the LMT training step, the simple logistic method is employed [34]. The aim of the simple logistic method is to control the parameter numbers, and to keep the model as simple and easy as possible. Therefore, new parameters are added, step by step, to improve the performance of each function at each node in the tree during training (see Section 3). This can also avoid the problem of model significance in the logistic regression, especially in the multiple logistic functions that build a full logistic model with all parameters. However, only some of the functions have less parameters than selected, indicating that most of our selected parameters affect predictive performance.         Many metrics can be employed to assess the predictive results of classification problems. We employed accuracy and Cohen's Kappa values for the evaluation herein. For a confusion matrix with n rows and columns, i is the row counter, j is the column counter, and m is the element in the matrix. The accuracy can be computed as [51]:

No. Stability Condition Regression Models
The Cohen's Kappa coefficient can also be employed to assess predictive performance. However, it can only be employed on binary prediction results, so only Model 2 is evaluated by this metric. It is a robust measure of the proportion of cases classified correctly after the probability of chance agreement has been considered and removed [52]. The expression is [53]: where Ac is the accuracy value obtained in Equation (5), and P e is the hypothetical probability of chance agreement, defined as Equation (7).
where m i+ and m + i are the sum of elements in the i-th row and i-th column of the confusion matrix. The results of κ reflect the agreement condition of the results, with a complete agreement corresponding to a value of 1 and no agreement corresponding to less than 1. The accuracy for Model 1 was equal to 94.1%, showing a satisfying predictive performance. For Model 2, accuracies for the training results with and without cost matrix were 92.9% and 93.1%, respectively, which conform to Kappa coefficients of 0.80 and 0.79. Both metrics indicate that the predictive results are reliable and acceptable.

Cross-Validation
A 10-fold cross-validation procedure was adopted to validate Models 1 and 2 and to further examine their predictive results. For 10-fold cross-validation, each database was randomly sorted into 10 datasets; then, for each dataset, the remaining nine datasets were used to train the model, and the initially selected set was used for pillar stability evaluation with the trained LR models and to examine their performance using the observations within the set. The process was then repeated for all the ten sets. Results are listed in Tables 11 and 12, revealing that the average accuracy values for Models 1 and 2 are 79.1% and 80.5%, respectively.

Validation through New Cases
Validation of model performance can also be conducted using new case histories. Therefore, another nine and twelve case histories (for Models 1 and 2, respectively) that were not employed for training work were used for validation. Their observed pillar stability conditions, as well as the predictive results obtained with functions in the tree models, are reported in Tables 13 and 14, where one can observe that our proposed tree models performed very well for most cases. For Model 1, all the case histories are well predicted. For Model 2, only three cases out of 12 are wrongly predicted. The test results confirm that the two proposed logistic tree models can reliably predict new case histories.

Comparison
To further assess our proposed Models 1 and 2, the results of several machine learning methods that have been studied in the literature, and of some other techniques that have not yet been employed for pillar stability prediction, are compared with the results of our models. Table 15 compares prediction results for Model 1. It should be noted that although all the models proposed in the literature achieved good predictions with different input parameters, our proposed model has superior predictive performance. We can also note that the random forests shows the optimal predictive performance among all techniques (Techniques No. 6 to No. 12 are trained by WEKA); however, its cross-validation result is lower than the model proposed in this research. In other words, our proposed Model 1 has better predictive performance in both the training and cross-validation stages.  Table 16. Note that fandom forests, decision trees and our proposed LMT show better predictive performance than the other methods considered; however, when compared to the other two algorithms (random forest and decision trees), our proposed model has the advantage that it can also be utilized for probabilistic risk analyses (discussed in the next section). In addition, the main advantage of the proposed Models 1 and 2 is that they illustrate an explicit class probability using input parameters, which can be easily employed to predict the stability and potential risk of pillars using simple calculation by engineers. (The developed MATLAB "m-files" are provided in the Supplementary Material File S1). However, in the "black box" nature of other machine learning techniques, the training and predicting procedures can only be accomplished using other tools such as WEKA, R, or MATLAB. Therefore, they cannot be easily applied for quick forecast directly. For instance, the random forests algorithm may predict pillar stability more accurately, but it does not provide any equations that can be used for such evaluation.

Risk Analysis
We use two case histories in Table 13; Table 14 to present an application of the proposed Models 1 and 2 for risk analyses. The group of typical input data for Case 1 is as follows: w = 3.5 m, h = 3.8 m, r = 0.92, UCS = 94 MPa, and p = 55 MPa; for Case 2, we have H = 106.68 m, w = 12.19 m, B = 6.1 m, and h = 4.27 m. Then, given the tree structures from Figures 3 and 4, we employ the proper functions to calculate the probability of failure (PoF). For instance, as the value of p for Case 1 is 55 MPa, we take the right branch of the tree for Model 1, and then we check the r value. With the r value being equal to 0.92, we turn to the left branch and compare the p value again; then the final function (LM2) for this case is decided by the UCS value of 94 MPa. According to the LM2 equations in Table 5, we first calculate function values (F S , F U , and F F ); then, using Equation (3), the probability of rock pillar stability can be expressed by Equation (11): Finally, we can obtain the results: Ps = 0, P U = 37.8%, and P F = 62.2%, which implies that this case has a PoF equals to 62.2%. For Case 2 of Model 2, the function (LM6) in Table 6 can be used to calculate the PoF using the same method. Results show that the probability of a stable coal pillar for Case 2 is 91.1%, which also agrees with the actual conditions. Such probability results can be incorporated into risk analyses with the corresponding failure cost assessment later.

Feature Importance Analyses
In the training process of Models 1 and 2, information gain measures the separation of the training cases according to their target classification by a given attribute. It is employed to select among the candidate features during the tree growing procedure. We suppose that a case history is (x, y) = (x 1 , x 2 , x 3 ,..., x k , y) and x k is the k-th attribute of that case with the corresponding class label y. The expected information gain is the change in information entropy H from an initial state to a state that takes some information, computed as [54]: where S denotes the training data set and S t is the subset of S for which feature k has value of t (S t = {x ∈ S|x k = t }). In the tree model training, before determining the feature in the root node, all the information gains provided by different features are calculated and compared, and the most influential feature is chosen as the root node in the tree structure. The other features in the nodes of the decision tree appear in descending order of importance [32,55]. Compared to other machine learning techniques, for which extra analyses are required to analyze sensitivity, the LMT decides the most influential feature during the tree growing process. Therefore, the average stress p is the most significant parameter in the pillar stability prediction for Model 1, followed by the ratio of pillar width to height r. Similarly, the most critical parameter for Model 2 is also r, followed by the excavation depth H. Both importance analyses show that pillar stability is largely influenced by the applied stress, ratio of width to height, and depth, thus agreeing with the empirical observations mentioned earlier in this paper.

Conclusions
A novel application of Logistic Model Trees (LMT) for instability assessment of underground pillars is presented. Tree structures and the corresponding functions are employed to assess the stability of pillars, given information on several features-including pillar width, pillar height, ratio of pillar width to height, uniaxial compressive strength of rock, average pillar stress, underground depth, and bord width-based on which predictions are conducted. The LMT is learned by LogitBoost, using two sets of case histories from the databases that we compiled from the literature. Results show that both Models 1 and 2 can accurately predict the stability of rock and coal pillars.
The trained tree models were validated using the original databases and 10-fold cross-validation. In addition, nine and twelve new cases-not previously employed for training-were used to further validate the models proposed. Moreover, predictions of several other machine learning methods were also compared with our proposed models. Results showed that the accuracies of the proposed tree models are among the highest and that they are probably adequate for mine site applications. In addition, results suggest that the LMT technique can offer useful information about the probability of pillar failure during underground excavation, so it can be utilized for risk analyses of pillar stability.
Furthermore, feature importance analyses were conducted to appraise the most influential input features for pillar stability. The average pillar stress p was found as the most influential parameter for Model 1, with other factors such as r and UCS also having a critical influence on predictions. For Model 2, r showed the highest effect on the output, followed by H and w. It is expected that these results can lead mine engineers to focus their efforts towards those crucial identified factors.
As the dataset employed for LMT learning and training was limited, it is expected that its predictive capacity could be increased if more case histories are collected. In addition, there are a few minor influencing parameters such as pore-water pressure, gas and water contents of the rocks, etc., which are neglected in this study due to the unavailability of data, and could have been among the selected parameters. More importantly, the application of the proposed models should be within the range of selected parameters.
Finally, the main advantage of LMT compared to other related algorithms commonly employed for pillar stability estimation is that it can be trained easily (even with more input parameters) and that its tree structure, with a LR function in each leaf, explicitly demonstrates the relationship between the inputs and predictive outputs. In addition, it is expected that, with its intuitive features and easy implementation, the LMT can also be employed to solve other geotechnical problems in the future.