Landslide Susceptibility Mapping Based on Random Forest and Boosted Regression Tree Models, and a Comparison of Their Performance

: This study aims to analyze and compare landslide susceptibility at Woomyeon Mountain, South Korea, based on the random forest (RF) model and the boosted regression tree (BRT) model. Through the construction of a landslide inventory map, 140 landslide locations were found. Among these, 42 (30%) were reserved to validate the model after 98 (70%) had been selected at random for model training. Fourteen landslide explanatory variables related to topography, hydrology, and forestry factors were considered and selected, based on the results of information gain for the modeling. The results were evaluated and compared using the receiver operating characteristic curve and statistical indices. The analysis showed that the RF model was better than the BRT model. The RF model yielded higher speciﬁcity, overall accuracy, and kappa index than the BRT model. In addition, the RF model, with a prediction rate of 0.865, performed slightly better than the BRT model, which had a prediction rate of 0.851. These results indicate that the landslide susceptibility maps (LSMs) produced in this study had good performance for predicting the spatial landslide distribution in the study area. These LSMs could be helpful for establishing mitigation strategies and for land use planning.


Introduction
A landslide is defined as a natural disaster that occurs when gravity causes a mass of debris, soil, or rock to move on a downward slope [1]. The majority of landslides occur as a result of hydroclimatic events, such as prolonged or intensive rain. Furthermore, mechanisms such as seismic triggers, wind, and freeze-thaw cycles are known to initiate landslides [2].
Mountains with shallow layers of soil that have formed in place from weathered gneiss and granite make up roughly 70% of the Korean peninsula [3]. Such terrain is vulnerable to weakening during heavy rainfall. Most of the annual precipitation occurs during the summer, when heavy rain and typhoons frequently occur. In particular, the heavy rain associated with typhoons has the potential to cause landslides in South Korea [4]. The year 2011 was a particularly devastating year, with 43 landslide-related casualties in Chuncheon and at Woomyeon Mountain in the area surrounding Seoul City. This is the largest number of landslide-related casualties since 2000.
South Korea has not been alone in experiencing an increase in such natural disasters. Other regions around the world have also experienced more frequent landslides on a larger scale and with more severe damage. In future decades, this trend will probably continue because of ongoing deforestation, increased urbanization, and an increase in regional precipitation in landslide-prone areas due to climate change [5]. It is essential that both susceptible and stable areas be identified to mitigate property damage, environmental degradation, and loss of life. Consequently, landslide susceptibility assessments, i.e., assessments of the spatial probability of a landslide occurring, are a huge step forward in the comprehensive hazard management of landslides [6,7]. The landslide susceptibility map (LSM) produced by a landslide susceptibility assessment can be a useful tool for authorities with decision-making capabilities.
Many methods and techniques have been proposed to evaluate landslide susceptibility. In the past few decades, statistical approaches have become popular in the use of remote sensing (RS) with a geographic information system (GIS). There are many statistical approaches used in landslide susceptibility assessment, including a frequency ratio (FR) [8,9], certainty factor (CF) [10], statistical index (SI) [11,12], as well as weight of evidence (WoE) [7,13,14] and logistic regression (LR) [15,16] approaches.
Recently, machine learning techniques have become popular in various fields. Machine learning, a branch of artificial intelligence, uses computer algorithms to analyze and predict information based on learning from training data [17,18]. Due to its robustness and high generalization capability, the use of machine learning has increased in landslide susceptibility analysis. Among the machine learning methods, artificial neural network [19,20], fuzzy logic [21,22], neuro-fuzzy [23], support vector machine [24,25], random forest [26,27], and naïve Bayes tree [17,28] methods have been popularly applied.
More recently, ensemble machine learning techniques have been used to enhance the prediction power and robustness of landslide susceptibility assessment. The ensemble methods, formed by a combination of variously based classifiers, have typically demonstrated significant improvement [17,24,29,30]. Ensemble techniques, which are relatively new approaches for producing a landslide susceptibility map, have been rarely used in the field. Therefore, the main objective of this research was to analyze and compare the performance of different ensemble models-namely, the random forest (RF) and boosted regression tree (BRT) models-for landslide susceptibility analysis. The RF and BRT models are very popular ensemble methods. Both are tree-based algorithms that predict the results by combining individual trees. However, the RF and BRT models build trees in different ways. Considering these characteristics, these models are appropriate for producing LSMs and for comparing LSM results. The results of the models were compared using the receiver operating characteristic (ROC) curve and statistical indices to determine the more robust model.

Study Area
The study area, Woomyeon Mountain, is located in the Seocho district of Seoul City, South Korea. This area lies within 37 • 27 00 -37 • 28 55 N and 126 • 59 02 -127 • 01 41 E (Figure 1). The average elevation is 293 m above sea level, and the slope is approximately 30 • -35 • . The bedrock is Precambrian banded biotite gneiss, which is believed to be highly susceptible to landslides because of severe weathering and abundant faults ( Figure A1). In addition, granite gneiss with relatively poor compositional differentiation is excavated en masse, and there is partial distribution of an embedded dike. The gneiss outcrop is poor, as a result of severe weathering in the overall area, and its foliation structure is irregular, due to several folding events [31].
This area experienced concentrated precipitation from 26-29 July 2011. The maximum precipitation, which occurred during 2 h one morning, was 164 mm. This exceeded the 156-mm, 100-year return period. This heavy precipitation led to a debris flow landslide in the area near Woomyeon Mountain, and 1-1.5 m of stratum flowed over areas near the mountain. Seven locations in the study area, including two locations in the valley area damaged primarily from flooding and five locations damaged by debris flow, were affected by the landslide. The total area damaged by the debris was approximately 276,683 m 2 , and the maximum length of damage from the upper part of the steep-slope disaster area to diffuse areas was approximately 764 m. The event caused 16 deaths and 10 building collapses [31].
Redlands, CA, United States). Among the digitized landslide locations, landslide locations belonging to rupture zones were converted to point data using a centroid technique. The point data representing the landslide locations were converted to a pixel format, with resolution of 10 m. From the 140 identified landslides, 42 (30%) were reserved to validate the model, after 98 (70%) had been chosen at random for model training. Additionally, non-landslide pixels were selected randomly from the non-landslide area: 98 non-landslide pixels were used for the training dataset, and 42 non-landslide pixels were used to build the validation dataset. This generating and splitting process was performed repeatedly more than 10 times. Finally, the combination utilized was found through the area under the receiver operating characteristic (ROC) curve (AUC) method.

Landslide Inventory
Landslide locations were identified using 32 aerial photographs of the study area, taken after the occurrence of the landslides. These aerial photographs were taken by a digital mapping camera with a spatial resolution of 10 cm. The orthorectified photographs were produced using the Leica Photogrammetry Suite (LPS) mounted on ERDAS Imagine 2011 (Erdas, Inc., Norcross, GA, United States). Landslide locations were digitized by visual interpretation using ArcGIS 10.2 (ESRI, Inc., Redlands, CA, USA). Among the digitized landslide locations, landslide locations belonging to rupture zones were converted to point data using a centroid technique. The point data representing the landslide locations were converted to a pixel format, with resolution of 10 m. From the 140 identified landslides, 42 (30%) were reserved to validate the model, after 98 (70%) had been chosen at random for model training. Additionally, non-landslide pixels were selected randomly from the non-landslide area: 98 non-landslide pixels were used for the training dataset, and 42 non-landslide pixels were used to build the validation dataset. This generating and splitting process was performed repeatedly Appl. Sci. 2019, 9, 942 4 of 19 more than 10 times. Finally, the combination utilized was found through the area under the receiver operating characteristic (ROC) curve (AUC) method.

Landslide Explanatory Variables
Landslides usually occur by complex interactions among various explanatory variables, and there is no consensus about which landslide explanatory variables to use. In this study, 14 explanatory variables were selected, based on a literature review and data availability. These factors were divided into the following three categories: topography, hydrology, and forestry (Table 1, Figure A2). These factors were produced in raster format with a cell size of 10 × 10 m, considering the scale of the input data, using ArcGIS 10.2 and ERDAS Imagine 2011; the total number of cells in the study area was 67,005. For the next process, the continuous variables among the explanatory variables were reclassified into seven classes, using ArcGIS 10.2. ArcGIS 10.2 provides various classification schemes, such as equal interval, standard deviation, natural break, quantile, etc. Natural break classification groups the classes based on break points that are relatively large jumps in data values. This classification method can be used to maximize the variance between classes. In addition, Cao et al. (2016) [32] indicated that natural break classification is more appropriate for the classification of variables, because their results showed that the LSM produced had higher accuracy compared to that using a different classification method. Therefore, natural break classification was used in this study. Topography factors include altitude, slope degree, slope aspect, profile curvature, and plan curvature. Altitude is an influential factor among the various landslide explanatory variables, because it is affected by several geomorphologic and geological processes. Slope, which can be described as the form between any section of the surface and a horizontal datum, has considerable influence on slope stability [33]. The degree of vulnerability to landslides may differ based on slope direction, because the water content of the surface, vegetation type, and soil strength may be different. In addition, both the profile and plan curvatures can be classified as flat, concave, or convex. During the rainy season, concave slopes may contain more moisture than convex slopes or flat slopes, so the concave slopes may be more vulnerable to landslides. All of these variables were extracted from the 10-m digital elevation model (DEM), using the spatial analyst tool of ArcGIS. The DEM was produced from 1:5000 topographic maps provided by the Korean National Geographic Information Institute.

Hydrology Factors
The hydrology factors were distance to streams, topographic wetness index (TWI), stream power index (SPI), sediment transport index (STI), and terrain roughness index (TRI). The streams were delineated by flow accumulation and converted to a vector format. The distance to streams was calculated using the Euclidean distance function in ArcGIS. Beven and Kirby (1979) [34] developed a TWI that reflects water's tendency to accumulate anywhere within the catchment area, accumulations that will then tend to move downslope as a result of gravity [35]. The water flow's power to erode is measured by the SPI, based on the assumption of proportionality of discharge to a catchment's specific area [36]. The STI is also often used to reflect the overland flow's power to erode [37]. The TWI, SPI, and STI were calculated with their base in specific catchment areas (A s ) and slope maps, using the following: where A s represents the specific catchment area (m 2 /m), and β represents the local slope gradient (degrees).
In addition, the TRI, which represents the concave and convex upward slopes [38], was calculated as where max and min represent the maximum and minimum values of altitude among the nine rectangular neighbor pixels, respectively.

Forestry Factors
Vegetation prevents erosion on a slope by buffering the impact of rain falling on the slope, and vegetation roots increase the shear strength of the slope by increasing the shear strength of the soil. The forestry factors include timber type, timber diameter, timber age, and timber density. Here, timber type and timber age mean the species and average age of planted trees, respectively. In addition, timber diameter represents the size of the diameter at chest height. Timber density refers to the degree of closure of the crown canopy. These values were obtained from a 1:5000 scale forest map produced by the Korea Forest Research Institute.

Methodology
This study was performed using the following main steps: (1) collection and construction of database of landslides and landslide explanatory variables, (2) preparation of the training and test datasets through repeated random sampling, (3) feature selection using information gain (IG), (4) landslide susceptibility mapping using RF and BRT models, and (5) validation and comparison of performance among landslide susceptibility maps (LSMs) (Figure 2). The IG, RF, and BRT models were implemented in R (Foundation for Statistical Computing, Vienna, Austria) using the "FSelector," "randomForest," and "gbm" packages, respectively. These algorithms were performed employing a 10-fold cross-validation approach, to reduce the variability of the model results. Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 18

Landslide Dataset Preparation
In this study, correlations between landslide and landslide explanatory variables were analyzed using the FR (Table A1). The FR is the ratio of the area where landslides occurred to the total study area. The FR was calculated by dividing the ratio of landslide occurrence, for the class or type of each factor, into an area ratio, for the class or type of each factor to the total area. The calculated FRs pertaining to each landslide explanatory variable were normalized from 0 to 1. The normalized FRs were extracted for the landslide and non-landslide dataset. Subsequently, these data were used for the training and validation datasets, to run the models and evaluate the prediction capabilities of the models, respectively.

Information Gain
The landslide explanatory variables have a crucial role in producing LSM. Some landslide explanatory variables might be associated with reductions in model performance, overfitting, model training time, and predictive capability [39]. Therefore, it is necessary to recognize and choose proper landslide explanatory variables.
Various methods such as IG [17], chi-square statistics [30], and Relief-F [29] have been proposed for feature selection in landslide modeling. In this study, the IG, proposed by Quinlan (1993) [40], was used to determine irrelevant and unimportant variables. The IG evaluates an attribute by determining the overall information gain in terms of the class. Consequently, the result can determine the ranking of importance, based on the normalized average merit contributed by each attribute [41,42].
The IG value of landslide explanatory variable belonging to class L (landslide and nonlandslide) is calculated as [24,30] where ( ) is the entropy value of L, and ( | ) is the entropy of L after integrating the values of landslide explanatory variable . These values are calculated using Equation (6) and Equation (7), respectively: where ( ) is the prior probability of the class L, and ( | ) is the posterior probabilities of L given the values of explanatory variable . An explanatory variable with a higher IG value has higher rank, meaning that it is more important to landslide models. By contrast, an explanatory

Landslide Dataset Preparation
In this study, correlations between landslide and landslide explanatory variables were analyzed using the FR (Table A1). The FR is the ratio of the area where landslides occurred to the total study area. The FR was calculated by dividing the ratio of landslide occurrence, for the class or type of each factor, into an area ratio, for the class or type of each factor to the total area. The calculated FRs pertaining to each landslide explanatory variable were normalized from 0 to 1. The normalized FRs were extracted for the landslide and non-landslide dataset. Subsequently, these data were used for the training and validation datasets, to run the models and evaluate the prediction capabilities of the models, respectively.

Information Gain
The landslide explanatory variables have a crucial role in producing LSM. Some landslide explanatory variables might be associated with reductions in model performance, overfitting, model training time, and predictive capability [39]. Therefore, it is necessary to recognize and choose proper landslide explanatory variables.
Various methods such as IG [17], chi-square statistics [30], and Relief-F [29] have been proposed for feature selection in landslide modeling. In this study, the IG, proposed by Quinlan (1993) [40], was used to determine irrelevant and unimportant variables. The IG evaluates an attribute by determining the overall information gain in terms of the class. Consequently, the result can determine the ranking of importance, based on the normalized average merit contributed by each attribute [41,42].
The IG value of landslide explanatory variable C i belonging to class L (landslide and non-landslide) is calculated as [24,30] where IF(L) is the entropy value of L, and IF(L|C i ) is the entropy of L after integrating the values of landslide explanatory variable C i . These values are calculated using Equation (6) and Equation (7), respectively: Appl. Sci. 2019, 9, 942 7 of 19 where P(L i ) is the prior probability of the class L, and P(L i |C i ) is the posterior probabilities of L given the values of explanatory variable C i . An explanatory variable with a higher IG value has higher rank, meaning that it is more important to landslide models. By contrast, an explanatory variable with an IG value of zero must be removed from the dataset, because that factor does not make a contribution.

Landslide Susceptibility Analysis
3.3.1. Random Forest RF, developed by Breiman (2001) [43], is a popular ensemble learning method that has been used widely for classification, regression, clustering, and interaction detection. A single decision tree is a weak classification, because of its high variance and bias. However, RF tends to produce robust models, because it can mitigate these problems by using ensemble trees [44].
RF generates thousands of random binary trees to form a forest. Each tree is grown based on a bootstrap sample, using a classification and regression trees (CART) procedure with a random subset of variables selected at each node [26,45]. For each tree grown on a bootstrap sample, the "out-of-bag" (OOB) error rate is calculated using observations left out of the bootstrap sample. The final decisions of class membership and model construction (output) are determined by the majority vote among all trees [46].
Two types of error rate-the mean decrease in accuracy and the mean decrease in the Gini coefficient-were calculated. These measures have been widely used to rank and select variables [26,47].
To run the RF model, the user should optimize two priori parameters, the number of trees in the forest (ntree) and the number of variables tested at each node (mtry), to minimize the OOB error and obtain good model performance [44,45].

Boosted Regression Tree Model
The BRT model is a combination of statistical and machine learning techniques. The BRT model fits different techniques and combines them to improve the performance of a single model [48,49]. Two different algorithms, namely boosting and regression, are used in the model, and the strengths of these algorithms are combined to improve model accuracy and decrease model variance [45,50]. Boosting is one of the most powerful learning methods for improving model accuracy, by iteratively fitting new trees to the residual errors (RE) of the existing tree assemblage [45,51]. In addition to boosting, the BRT model uses regression trees in the modeling process. Regression trees are categorized from the classification and regression tree approaches from the decision tree group of models [52].
In the model, among the various parameters, the number of trees is automatically set through internal cross-validation. In addition, the learning rate, the number of nodes in a single tree, and bag fraction were determined through a trial-and-error approach [53]. The complexity of the model and the contribution of each tree to the model are controlled by a shrinkage parameter and the learning rate, respectively. The bag fraction and shrinkage parameter determine the number of trees required to reach the optimal solution [54].

Confusion Matrix
The confusion matrix includes true positive (TP), false positive (FP), true negative (TN), and false negative (FN) categories. Using these values, various statistical indices, such as accuracy, sensitivity, specificity, threat score, equitable threat score, Pierce's skill score, odds ratio, and odds ratio skill score can be calculated [55]. The value calculated from the confusion matrix provides useful information on model performance and classification accuracy.
In this study, the sensitivity, specificity, overall accuracy, and kappa statistic were used to validate the performance of the LSMs. The percentages of landslide and non-landslide pixels classified correctly into those two categories enable the calculation of sensitivity and specificity, and the overall percentage classified correctly (in both categories together) indicates the accuracy of the LSMs [56]. In addition, the kappa statistic is used to evaluate the reliability of the landslide models. Its value ranges from −1 (non-reliable) to 1 (reliable) [57].

Receiver Operating Characteristic
The receiver operating characteristic (ROC) curve has been commonly used to validate the quality of a probabilistic model. The ROC curve is plotted by statistical index value pairs, with the false positive rate (sensitivity) on the x-axis and the "100−false negative rate" (100−specificity) on the y-axis. The ROC curve can be classified as a success rate curve or prediction rate curve, depending on the dataset used. The success rate curve, calculated using the training dataset, represents how well the LSMs fit the data. The prediction rate curve, calculated using the validation dataset, represents how well the model and landslide explanatory variables predict a landslide [11]. The ROC curve can be verified quantitatively when the area under the ROC curve (AUC) is calculated. AUC values range from 0.5 to 1.0. AUC values closer to 1 indicate a more accurate model.

Selection of Landslide Explanatory Variables
The average information gain (AIG) value, and its standard deviation for each landslide explanatory variable, were calculated and ranked ( Table 2). All landslide explanatory variables used in this study contributed to the landslide models, because the AIG values of these variables were more than 0. According to the results, the TRI had the highest AIG value (0.086), which means that this factor made the greatest contribution to the landslide models in this study area. By contrast, timber diameter made the smallest contribution to the landslide models, as indicated by the lowest AIG value (0.005).

Training the Random Forest and Boosted Regression Tree Models
The training dataset was used to train the RF and BRT models for landslide susceptibility assessment. During the training process, the optimum values of the parameters for the models were applied to obtain high model predictive capability. The optimized values for the RF model were 300 for ntree and 2 for mtry. In the case of the BRT model, the optimized values for n.trees, interaction.depth, shrinkage, and n.minobsinnode were 500, 1, 0.01, and 10, respectively. Subsequently, the RF and BRT models were constructed using the optimized parameters, based on the training dataset.
After their construction, the RF and BRT models were applied throughout the whole study area to produce LSMs.

Model Validation and Comparison
The performance of each model was analyzed using the training dataset. The RF model showed a higher sensitivity value (98.00%) than did the BRT model (79.57%). This result showed that the RF model classified more correctly than the BRT model in the landslide class. The specificity results also indicated that the RF model had higher specificity (100.00%) in the non-landslide class, indicating that the non-landslide pixels were more correctly classified. The specificity value of the BRT model was 76.70%. Because of the lower sensitivity and specificity values of the BRT model, the overall accuracy and kappa index values were lower, with values of 78.16% and 0.561, respectively. In the case of the RF model, the overall accuracy and kappa index were 98.98% and 0.980, respectively. In addition, the success rate and the prediction rate were analyzed using the training dataset and the validation dataset, respectively (Figure 3). In the case of success rate, the RF and BRT models had values of 0.999 and 0.887, respectively. The prediction rate curve also showed that the RF model had a higher AUC (0.865) than the BRT model (0.851). Overall, the AUC values of all models were greater than 0.8. These results show that the LSMs constructed in this study have good accuracy in the spatial prediction of landslide susceptibility.

Model Validation and Comparison
The performance of each model was analyzed using the training dataset. The RF model showed a higher sensitivity value (98.00%) than did the BRT model (79.57%). This result showed that the RF model classified more correctly than the BRT model in the landslide class. The specificity results also indicated that the RF model had higher specificity (100.00%) in the non-landslide class, indicating that the non-landslide pixels were more correctly classified. The specificity value of the BRT model was 76.70%. Because of the lower sensitivity and specificity values of the BRT model, the overall accuracy and kappa index values were lower, with values of 78.16% and 0.561, respectively. In the case of the RF model, the overall accuracy and kappa index were 98.98% and 0.980, respectively. In addition, the success rate and the prediction rate were analyzed using the training dataset and the validation dataset, respectively ( Figure 3). In the case of success rate, the RF and BRT models had values of 0.999 and 0.887, respectively. The prediction rate curve also showed that the RF model had a higher AUC (0.865) than the BRT model (0.851). Overall, the AUC values of all models were greater than 0.8. These results show that the LSMs constructed in this study have good accuracy in the spatial prediction of landslide susceptibility.
(a) (b) Figure 3. Analysis of the receiver operating characteristic (ROC) curve for the two landslide susceptibility maps: (a) success rate curve using the training dataset and (b) prediction rate curve using the validation dataset.

Generating Landslide Susceptibility Maps
The RF and BRT models were used to develop LSMs in the study area. The LSMs were prepared by generating landslide susceptibility indices (LSIs) and reclassifying the class. The LSIs were calculated based on the trained RF and BRT models. Using the natural breaks method, the LSMs were reclassified into five susceptibility classes: very high, high, moderate, low, and very low (Figure 4). Overall, the distribution of LSI for each susceptibility class was similar between the LSM produced by RF (RF LSM) and that produced by BRT (BRT LSM). The "high" and "very high" susceptibility classes covered about 30% of the total area. The RF model had a value of 34.69%, and the BRT model had the lower value of 31.11%.

Generating Landslide Susceptibility Maps
The RF and BRT models were used to develop LSMs in the study area. The LSMs were prepared by generating landslide susceptibility indices (LSIs) and reclassifying the class. The LSIs were calculated based on the trained RF and BRT models. Using the natural breaks method, the LSMs were reclassified into five susceptibility classes: very high, high, moderate, low, and very low (Figure 4). Overall, the distribution of LSI for each susceptibility class was similar between the LSM produced by RF (RF LSM) and that produced by BRT (BRT LSM). The "high" and "very high" susceptibility classes covered about 30% of the total area. The RF model had a value of 34  The LSMs produced from the two models were validated based on the landslide density (LD) of each susceptibility class on the LSMs. The LD is the ratio of the percentage of landslide pixels to the percentage of all pixels for each susceptibility class shown on the map [56]. LD was calculated by overlaying the five LSMs and the landslide inventory map. Generally, for the study area, the value of LD increased gradually, from very low to very high susceptibility (Table 3). At the "very high" class, the RF and BRT models had LD values of 3.799 and 2.721, respectively. Overall, the models used in this study are suitable for LSM. Table 3. Landslide density on landslide susceptibility maps produced from the different models.

Pixels of Landslide
Landslide Density

Discussion
The LSMs produced using the models were evaluated by statistical indices and ROC curves. The RF model had better sensitivity, specificity, overall accuracy, and kappa values. The AUC values of the LSMs used in this study were about 80%, indicating reasonable accuracy. The RF model had higher AUC values for the success rate and prediction rate curves than the BRT model. Thus, these models had very high predictive performance. Furthermore, the LSMs would be produced differently depending on the methods used and the landslide explanatory variables selected. The landslide explanatory variables may not make equal contributions, which can affect prediction ability. In this study, the landslide explanatory variables used made different contributions to the models. Table 4 illustrates the importance of each explanatory variable, calculated and normalized in the RF and BRT models. In general, TRI had the highest importance to the models, whereas timber diameter, timber age, and timber density had lower predictive capability. The LSMs produced from the two models were validated based on the landslide density (LD) of each susceptibility class on the LSMs. The LD is the ratio of the percentage of landslide pixels to the percentage of all pixels for each susceptibility class shown on the map [56]. LD was calculated by overlaying the five LSMs and the landslide inventory map. Generally, for the study area, the value of LD increased gradually, from very low to very high susceptibility (Table 3). At the "very high" class, the RF and BRT models had LD values of 3.799 and 2.721, respectively. Overall, the models used in this study are suitable for LSM. Table 3. Landslide density on landslide susceptibility maps produced from the different models.

Random Forest
Boosted Regression Tree

Pixels of Landslide
Landslide Density

Discussion
The LSMs produced using the models were evaluated by statistical indices and ROC curves. The RF model had better sensitivity, specificity, overall accuracy, and kappa values. The AUC values of the LSMs used in this study were about 80%, indicating reasonable accuracy. The RF model had higher AUC values for the success rate and prediction rate curves than the BRT model. Thus, these models had very high predictive performance. Furthermore, the LSMs would be produced differently depending on the methods used and the landslide explanatory variables selected. The landslide explanatory variables may not make equal contributions, which can affect prediction ability. In this study, the landslide explanatory variables used made different contributions to the models. Table 4 illustrates the importance of each explanatory variable, calculated and normalized in the RF and BRT models. In general, TRI had the highest importance to the models, whereas timber diameter, timber age, and timber density had lower predictive capability. From the results, ensemble classification, such as that done by the model used in this study, can improve the performance of single (weak) classifiers and the prediction accuracy of LSM [56]. However, the models had overfitting problems, as indicated by the AUC values calculated using the training and validation datasets. The AUC values of the success rate curve were very high, almost reaching a value of 1, but the AUC values of the prediction rate curve were lower. Especially in the case of the RF model, the AUC value of the prediction rate was decreased by about 20%. This result showed that the RF model was trained excessively by the training data. This can be associated with poor generalization from training data and increased error for real data. Overfitting is a common problem affecting researchers performing machine learning and data mining. There can be many reasons for overfitting. However, in this study, the landslide explanatory variables used still included noise, despite the feature selection process. In addition, because the landslide area is very small compared to the non-landslide area, the model could not learn and predict the non-landslide area.

Conclusions
This study compared and analyzed landslide susceptibility at Woomyeon Mountain using different models. For this purpose, landslide-related spatial data consisting of a landslide inventory, and landslide explanatory variables were collected and prepared. The landslide inventory map was built using aerial photographs. The 14 landslide explanatory variables were constructed from spatial data collected by government organizations. These factors included altitude, slope degree, slope aspect, profile curvature, plan curvature, distance to streams, TWI, SPI, STI, TRI, timber type, timber diameter, timber age, and timber density.
The contribution of each landslide explanatory variable was evaluated using the average IG value with a 10-fold cross-validation approach. All of the landslide explanatory variables contributed to the models, because the IG values of all factors were greater than zero. Therefore, the landslide susceptibility analysis and mapping were performed with all landslide explanatory variables using the RF and BRT models. The RF and BRT models were implemented in R. A popular open-source software, R is helpful for statistical computing and data visualization [58]. The models were constructed using optimized parameters, and LSI was predicted over the study area.
The LSMs produced in this study may prove useful for decision makers, planners, and engineers in disaster planning to minimize economic losses and casualties. In a future study, the accuracy of the LSMs of this study could be enhanced by selecting more optimal landslide explanatory variables and solving the problem of overfitting.