Introduction

Landslides are a complex hazard that occur on highlands or mountain slopes, conditioned by several topographical aspects related to soil properties, geological structure, lithology and climatic conditions, slope morphology, soil cover, and water flow (Margottini et al., 2013). The landslide triggers are commonly rainfall (Y. Liu, Xu, et al., 2021), earthquakes (Pang et al., 2022), deforestation (García-Ruiz et al., 2017), and human activity that changes the effect of topography (Li et al., 2020).

Landslides cause serious impacts that threaten humans (Hakim et al., 2022; Panahi et al., 2020), and damages to natural resources. Thus, landslide hazard assessment has become a task of interest for decision-making by government entities and municipal and/or urban planning departments.

In general, in regions where urban developments, residential areas, and service infrastructure coincide with mountainous terrain, the risk tends to be high for the population and the economic costs may include relocation of communities, reconstruction of structures, and restoration of the quality of water sources. In many developing countries, where land occupation has generally been carried out without adequate planning and in a disorderly manner, the growth of urban areas occurs in landslide-prone zones.

Susceptibility mapping for landslide prediction is a GIS-based method involving correlation of previous landslides with possible driving factors to identify areas at risk of landslides (Hakim et al., 2022; Hakim & Lee, 2020). In recent years, studies of landslide susceptibility mapping have employed various probabilistic and statistical methods (A. Saha & Saha, 2020; Silalahi et al., 2019). Given the complexity of landslide prediction, many researchers have turned their attention to using hybrid ensemble approaches that combine machine learning methods with metaheuristic algorithms (Jaafari et al., 2019) or ensemble learning techniques (Bui et al., 2019; B. T. Pham et al., 2019).

Currently, landslide susceptibility maps allow the identification of areas prone to the occurrence of a mass removal event where the potential damage to people and infrastructure must be reduced or controlled, however, the accuracy of these approaches varies according to the quality of the data, the model approaches used, and the landslide inventories.

The purpose of this study is to articulate the statistical approach of the landslide conditioning factors, Machine Learning Algorithms (MLA) and GIS, evaluating a flexible and agile methodology to estimate the landslide susceptibility defining areas prone to the landslide occurrence incorporating the interpretability criteria by mean of SHAP values approach. The MLA were validated in a case study in the “La Liboriana” River basin, located in the Municipality of Salgar (Antioquia) in the Colombian mountains Andes where Landslide Susceptibility Maps (LSMs) were obtained. The results can be used for mapping regional landslides to develop strategies to minimize the loss of human lives, infrastructure, and natural environment.

Antecedents

Due to the high levels of landslides, a great dynamic has been generated worldwide in the study of the associated phenomena in an attempt to understand the physical and economic aspects related to mass movements (Hidalgo & Vega, 2021). The United Nations Office for Disaster Risk Reduction (CRED and UNDRR, 2021) reports that in 2020 the average annual economic losses were higher than those of the last two decades, which amount to US$ 151.6 billion, in addition to an increase in the phenomena triggered by weather conditions.

In comparison to the previous two decades, 2020 was higher than the annual average in terms of number of recorded events and the annual average of economic losses, which is US$ 151.6 billion. There were considerably fewer deaths compared to the annual average of 61,709 and fewer people directly affected compared to the annual average of 201.3 million people. However, in 2020 there were 26% more storms than the annual average of 102 events, 23% more floods than the annual average of 163 events, and 18% more flood deaths than the annual average of 5,233 deaths (CRED and UNDRR, 2021).

Landslides cause loss of human life every year. Laccase and Nadim (2009) report that at least 17% of all natural hazard deaths worldwide are caused by landslides. Human losses derived by landslides occur predominantly in developing countries. In contrast, developed countries such as the United States and Japan report few human losses, but high annual economic losses, estimated between 1 and 6 billion dollars (Ospina-Gutiérrez & Aristizábal-Giraldo, 2021). It is estimated that the direct and indirect costs of mass movements can be significant in terms of gross domestic product (GDP), even in developed countries (Figure 1).

Figure 1.

Impacts of landslides in terms of gross domestic product (GDP).

Source. Hidalgo and Vega (2021).

In the case of Colombia, in the period from 2006 to 2014, 21,594 emergencies due to natural events were reported in the country, an average of 2,399 events per year. Of these, 14,641 (67.8%) were concentrated in the period from 2011 to 2014. As a result of the events that occurred, 3,181 people were reported dead in Antioquia, with 586 deaths (414 of them due to floods and landslides).

According to data from the Information System of Mass Movements (SIMMA, in Spanish) of the Colombian Geological Service (SGC, in Spanish), in the period between 1900 and 2018 at least 30,730 landslides have occurred in Colombia, which have left a balance of 31,198 fatalities and economic losses of USD$ 654 million (Ospina-Gutiérrez & Aristizábal-Giraldo, 2021). According to the DESINVENTAR database, 10,438 landslides have been recorded in Colombia between 1921 and 2020, leaving almost 7,313 dead and disastrous results for the country's economic system.

Figure 2 shows a summary of records of the DESINVENTAR database up to 2019 in the Department of Antioquia. According to the data reported, landslides occur more frequently between the months of May to July and September to November, which coincide with the hydrogeological characterization of the study area. In total, 1,566 events were reported, which show that the areas with the highest occurrence are located in the southwest (462 records) and east (358 records) zones excluding the Aburrá Valley.

Figure 2.

Landslides records for the Department of Antioquia.

Landslide Susceptibility Assessment (LSA)

Landslide susceptibility is the landslide probability of occurrence in a specific area, based on local terrain conditions by interaction of landslide conditioning factors (LCF). Usually, information about landslide magnitude is not available. LSA permits to identify areas potentially affected without considering the time lapse which a landslide might occur, or its magnitude. Commonly, LSA is based on statistical relationships between past landslides and LCF. Future landslides will occur under the conditions that led to past landslides.

In recent years, GIS and remote sensing data have been used to conduct many studies of disasters in mountain regions. Several researchers have built their methodology analyzing data of past landslides, and tested it through unknown landslides events. Different methodologies have been applied to spatially assess landslide susceptibility assessment (LSA). The main methods can be divided into qualitative (known as knowledge-driven or heuristic) or quantitative (data-driven and physically-based) methods (Lima et al., 2022). Their applicability and limitations can be found in literature. The first methods are based on the expert judgment, and usually involve qualitative terms to represent the susceptibility zoning. The quantitative methods establish numerical relationships between LCF and landslide occurrence (Marjanović et al., 2019).

The knowledge-driven or heuristic approach considers a direct mapping methodology establishing a direct relationship between the occurrence of landslides and the LCF using a landslide inventory at regional scale. This category may include the subjective geomorphological method, Analytic Hierarchy Process (AHP), Fuzzy Logic, Weighted Overlay, among others (Ali et al., 2021; Kaur et al., 2023; Q. B. Pham et al., 2021; Sahana & Sajjad, 2017; Sur & Singh, 2019).

The deterministic approach is based on slope stability methods. It is generally only applicable in relatively homogeneous terrain conditions throughout the study area and the types of landslides are known. It requires a high degree of simplification of the intrinsic variables to be used at local scale. This category includes geotechnical methods such as the Newmark’s Method (Infinite Slope), Bishop's Method, Morgenstern Price Method, among others. Various studies have used deterministic models, as the Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability (TRIGRS) model (Ma et al., 2021; Marin et al., 2021), Shallow Slope Stability (SHALSTAB) model (Pradhan & Kim, 2015), Stability Index MAPping (SINMAP) model (Michel et al., 2014), and Steady-State Infinite Slope Method (SSIS) (Si et al., 2020).

Finally, the data-driven or statistical approach is an indirect susceptibility methodology. It involves statistical analysis of the combinations of variables that led to landslide occurrence in the past. All possible intrinsic variables or LCF are entered and crossed into a GIS for analysis with a landslide inventory in a bivariate or multivariate way at regional scale (Dahal et al., 2012). This category also includes evaluation methods such as Frequency Ratio, Evidence of Weights, Linear and Logistic Regressions, among others. These statistical models have their own advantages and disadvantages and they have been widely used for LSA, but no agreement has reached to select the best method for landslide susceptibility analysis (Sahana & Sajjad, 2017; Zhang et al., 2020). Machine Learning Algorithms (MLA) are included in this category, but this is yet a topic of debate (Merghadi et al., 2020).

Machine Learning (ML) corresponds to a subset of a large discipline called artificial intelligence, which seeks to emulate human behavior through computer algorithms. ML uses statistical methods to train machines from data, that is, from experiences. Specifically on the topic of LSA using MLA, progress has been made since the early 2000s (Merghadi et al., 2020). MLA have been applied to solve geotechnical engineering problems widely in recent years, mainly because MLA are based on historical data, and they are more objective than the expert system methods. Moreover, they do not need more detailed mechanical parameters compared with the mechanical model methods (Y. Liu et al., 2019).

The MLA with the earliest development corresponds to Logistic Regression (LR), Artificial Neural Networks (ANN), Support Vector Machine (SVM), and Decision Trees (DT) (Bragagnolo et al., 2020; Dou et al., 2019; Nhu et al., 2020; Sahin, 2020; Sun et al., 2020; Wang et al., 2015). Moreover, in recent years, development has focused on bagging methods, including the very popular Random Forest (RF), boosting methods such as AdaBoost and XGBoost (Bui et al., 2019; Chen et al., 2020; Huang et al., 2020; Q. B. Pham et al., 2021; Sahin, 2020).

MLA have been increasingly used in LSA as a result of the fact that they can learn the association between landslide occurrences and LCF without the requirements and assumptions for a statistical model. It has been observed that the predictive power of conventional statistical methods is relatively low. The conventional statistical methods cannot accurately analyze the complex interrelationships between different causative factors. Due to accuracy and very high predictive capability, MLA are getting more importance and attention in spatial analysis of landslides.

With the development of ensemble learning, bagging, and boosting methods are increasingly being used for classification and regression (Chen et al., 2020). In this study, seven supervised MLA were employed. Two regression algorithms (Logistic Regressions) and five decision tree algorithms: RPART, CTREE, RF, Ranger, XGBoost, briefly described below.

The binary logistic regression algorithm seeks to study the relationship between a dichotomous response variable (i.e., it takes only two possible outcomes presence/absence) and one or more explanatory variables, which can be of both qualitative and/or quantitative nature (Hosmer et al., 2013). Two logistic regression algorithms were used, the conventional stepwise selection method and the other where the Least Absolute Shrinkage Selector Operator (LASSO) method is used which imposes a penalty on the regression coefficients and selects variables (Tibshirani, 1996).

RPART is a supervised learning algorithm where a classification tree is obtained if the response variable is dichotomous (discrete) or a regression tree if the response variable is continuous. Initially this algorithm finds the explanatory variable that best divides the data into groups based on a rule. Then, for each of the partitions, the process is repeated. This process is done recursively until it is impossible to find a better partition. A relevant feature of this algorithm is that a variable used to separate the data is not used afterwards (Therneau & Atkinson, 1997). Since the resulting tree is very large and becomes tedious to interpret, pruning techniques are used to reduce its size (Strobl et al., 2009).

CTREE are a special type of decision trees where the choice of variables is made by assessing whether there is an association between the response variable and each of the explanatory variables. If the null hypothesis of independence is not rejected for the whole set of hypotheses, the recursive process is stopped. Otherwise, the level of association of each set of significant tests is quantified, allowing new splits in the tree to be generated sequentially (Hothorn et al., 2006). One difference of this algorithm with the other decision tree algorithms is that no pruning is done for its statistical support.

Random Forest (RF) is an ensemble of decision trees, which are then combined into a single robust model (Breiman, 2001). RF uses a technique called bagging to build and train the ensemble of decision trees, allowing to reduce variance problems, prediction bias and overfitting when working with large amounts of data. One of the advantages of this classification algorithm is that it can handle many input variables and identify the most significant ones (dimensionality reduction) (Liaw & Wiener, 2002).

The Ranger algorithm is a fast implementation of the RF algorithm, particularly for handling high dimensional data (Wright & Ziegler, 2017). The XGBoost is a supervised MLA (Chen & Guestrin, 2016). The main idea of this algorithm is to generate multiple decision trees sequentially (boosting) where each one takes the results of the previous one and thus to generate an increasingly robust model with better predictive power. This process is repeated until the best possible model is obtained (Y. C. Chang et al., 2018).

Some of the aforementioned MLA have an internal structure that causes difficulty in explaining and interpreting their results, except those of a linear nature known as glass-box models. To address and solve the black-box issue of some MLA, as tree-based ensembles, kernel-based models, and neural networks, at both global and local levels, it is necessary to build eXplainable Artificial Intelligence (XAI) based solution. This approach provides an identification of the LCF influencing mainly effective classification of landslide susceptibility. Interpretable MLA can overcome the limitations of complex MLA in interpreting landslide susceptibility, and in fact, currently, few studies use the SHAP method to interpret the susceptibility of rainfall-induced shallow landslides (Zhou et al., 2022).

Methodology

Study area

Figure 3 shows the area corresponding to the case study in the “La Liboriana” River basin, located in the municipality of Salgar, in the southwest of the Department of Antioquia, western branch of the Colombian Andes. This basin joins El Barroso River basin, and both drain water into Cauca River, one of the most important rivers in the country (Hidalgo & Vega, 2021).

Figure 3.

Location of the study area.

The study basin presents geomorphological, geological, and weather conditions that make it particularly susceptible to landslides and flash floods. The area with slope gradients exceeding 30° accounts for 67% of the total area. It has a humid tropical climate with a mean annual temperature of 22 °C. The rainfall regime is dominated by high interannual and intraseasonal variability with a mean annual rainfall of 3,073 mm; and monthly rainfall distributions show evident seasonal patterns with two rainy seasons (peaks in May and October) (Ruiz Vásquez & Aristizábal, 2018).

The geomorphology of the basin exhibits in the upper part a mountain region with a rugged morphology and very steep forested hillslopes. In the middle and lower zones, grasslands (pastures) and coffee plantations have already substituted forest. The basin exhibits grazing areas and urban development near the riverbanks. Geologically, it is composed predominantly by a Cretaceous sedimentary rock formation and an intrusive Miocene body. These rocks have been severely weathered in situ under the humid tropical climate forming well-developed saprolite and residual soils (Ruiz Vásquez & Aristizábal, 2018).

LSA using machine learning algorithms

Seven supervised MLA were employed, two regression algorithms (Logistic) and five decision tree algorithms (Recursive Partitioning and Regression Trees [RPART], Conditional Inference Trees [CTREE], Random Forest [RF], Ranger, and Extreme Gradient Boosting Algorithm [XGBoost]). The LSMs were produced for each MLA. The LSA study was developed considering five phases: data processing, feature selection, dataset splitting and resampling, evaluation of the susceptibility, and comparison of the performance of the algorithms, as shown in Figure 4.

Figure 4.

Adopted methodology for landslide susceptibility assessment (LSA).

Data processing and landslide conditioning factors.

Several Landslide Conditioning Factors (LCF) have been employed in literature to produce the landslide susceptibility maps. Slope, aspect, lithology, plan curvature, and drainage density are most extensively used. Initially, data was processed in SAGA GIS 7.7 ( http://www.saga-gis.org) software and then collected on 16 LCF which can be grouped into categories as terrain, geological, hydrological, and coverage factors. The factors considered in the study are shown in Table 1 and then briefly described. All LCF layers were stacked into a geographic database in raster format, using ArcGIS 10.8 version ( https://desktop.arcgis.com).

Table 1.

Dataset of Landslide Conditioning Factors.

Additionally, a summary of the descriptive statistics for the chosen LCF is shown in Table 2. For the qualitative LCF, the absolute frequency and percentage of each of the categories are presented while for the quantitative, six descriptive statistics are presented, minimum (Min), quartile 1 (Q1), median (Median), mean (Mean), quartile 3 (Q3), and maximum (Max).

Table 2.

Statistical Description of Most Significant LCF.

Frequency Ratio (FR) which considers the effects of conditioning factors on landslide occurrence (Z. Chang et al., 2020) was determined. FR values were calculated as the ratio between the percentage of landslides and the percentage of class of each LCF. An FR value greater than 1.0, indicates a higher correlation between landslide and conditioning factors; whereas an FR value that is lower than 1.0, suggests a lower effect on landslide.

A landslide inventory documented in Ruiz Vásquez and Aristizábal (2018) for the study basin was used. This landslide inventory was obtained from a multi-temporal analysis of satellite images and aerial photographs, showing a landslide area covering approximately 0.6 km² corresponding to 1% of the basin.

Feature selection.

Statistical methods of correlation and multicollinearity analysis were used. The variance inflation index (VIF) and the tolerance (TOL) were analyzed. VIF focuses on the standard error variations of LCF, which implies the lower the standard error, the lower the multicollinearity risk, which in turn means the increased likelihood that robust results should be obtained using this LCF (Merghadi et al., 2020). LCF with a VIF value less than 5 were used in modeling. LCF with high correlation or multicollinearity were excluded from the training phase of MLA, as they may generate noise in the modeling through erroneous system analysis.

For LCF selection according to relative importance, a Relief F Method provided in WEKA open source ( https://www.cs.waikato.ac.nz/ml/weka) was used. This method calculates a weight value (average merit) for each variable to quantify its relevance (Dang et al., 2020). Features with the weights exceeding a certain threshold are selected for analysis. Factors assigned zero weight have no contribution to landslide occurrence and therefore, must be removed from further analysis (Hong et al., 2018). This study applied a cross-validation strategy, a standard random 5-fold cross-validation for Relief-F attribute selection method.

Split and resampling dataset.

For the purposes of the training and validation stages of the MLA, a usual 70/30 split was adopted (Achour & Pourghasemi, 2020; S. Saha et al., 2021). In order to balance the input dataset, over sampling techniques can be used to increase the minority class, that is, the landslide class, and for this purpose, the method SMOTE (Synthetic Minority Over sampling TEchnique) (Qing et al., 2020) was implemented under the R Studio environment ( http://www.rstudio.com). This method is an oversampling technique that creates synthetic minority class data points to balance the dataset using a K-nearest neighbor algorithm.

Selection of the best performance MLA.

An alternative method for evaluating and determining the statistical significance of systematic pairwise differences between the MLA is the Wilcoxon signed-rank test. This non-parametric test, as described by Dou et al. (2019) and Merghadi et al. (2020), was employed with a significance level of α = 5%. The test is based on a null hypothesis that assumes an equality between models only if the p-value is greater than 0.05. Rejecting the null hypothesis indicates statistically significant difference in performance between a pair of models, thereby establishing its reliability.

Once it has been identified that the models used showed statistically significant differences, seven statistical metrics were used to measure performance and validate the trained MLA. These metrics are widely used in the machine learning literature and are briefly described below. Five metrics were calculated based on the binary confusion matrix (Agarwal, 2020; Ahamad et al., 2020). To do this, first the following measures were calculated, True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN). Then, the following performance metrics were calculated:

(1)

(2)

(3)

(4)

(5)

Another statistical metric used was Area Under Curve (AUC) of the Receiver Operating Characteristic (ROC). The AUC ranges from 0 to 1 and represents the predictive performance of the model (Fawcett, 2006). There are five categories to understand this value with respect to the level of accuracy: excellent (0.9–1.0), good (0.8–0.9), fair (0.7–0.8), poor (0.6–0.7), and fail (0.5–0.6) (Das & Lepcha, 2019; Rasyid et al., 2016). In addition, the probability-based log loss classification metric was also calculated and used to compare the performance of MLA.

These metrics were calculated for each of the fitted algorithms. The choice of the best algorithm to explain the occurrence of a possible landslide was carried out by comparing and counting the number of metrics in favor of each algorithm. The MLA used in this study were implemented using R software (R Core Team, 2022). They are all part of the “alookr” package. For description of this package see  https://CRAN.R-project.org/package=alookr.

Interpretation of MLA output with SHAP method.

In this work, the SHapley Additive exPlanations (SHAP) approach was adopted to understand the MLA output. This method is a game theoretic approach to explain the output of any MLA. It provides a unified framework to interpret predictions through calculating the Shapley values that provide the coherence of the explanations (Inan & Rahman, 2023; Kavzoglu et al., 2021). The SHAP method allows local interpretation using the Shapley value for each feature in a single sample to show the contribution of each feature to the predicted value. The SHAP method can be consistent by aggregating local explanations into global explanations and by separating single-factor effects from interaction effects (Zhou et al., 2022).

Results and Discussion

Statistical description of most significant LCF

A multicollinearity analysis was performed. Then, the Relief-F attribute selection method for most significance LCF was used. LCF with average merit (AM) values less than or equal to zero were not considered in the MLA training process since they were not statistically significant in explaining the response variable. Figure 5 shows the significant LCF and their respective AM, indicating their relative importance for explaining the landslide occurrence. These results were obtained using a standard random 5-fold cross-validation with Relief-F method. All these values are positive, and indicate that NDVI has the highest AM, followed by the factors Elevation, Curvature, and TWI. Landcover and Soil Type have the lowest AM values. Table 2 shows the values of the tolerance and variance inflation factor for multicollinearity diagnosis. All the conditioning factors previously selected by the Relief-F method have a tolerance greater than 0.2 and a VIF less than 5, indicating that there are no multicollinearity problems among the LCF.

Figure 5.

Average merit of landslide conditioning factors considered.

Spatial distribution of most significant LCF and frequency ratio

Figure 6 shows the spatial distribution of the chosen LCF for modeling and its value ranges. In general terms it can be noticed that the spatial distribution of the elevation factor shows the altimetric variation of the study area, showing a definite pattern with a marked trend of increasing low values from the east to higher values in the west. In addition, it is observed that most of the landslides have occurred at higher elevations. In fact, it can be noticed that the greatest occurrence of landslides in the basin corresponds to the elevation range between 3,000 and 3,500 m.a.s.l., where there are landforms with gradients of slopes greater than 45°.

Figure 6.

Spatial distribution of landslide conditioning factors used in modeling.

The curvature convex and straight (plan) present the greater frequency ratio (FR), due to this formation is favorable for landslide initiation according to used landslide inventory. High values of Terrain Ruggedness Index (TRI) were present in the landslide zone where predominantly silty clay soils exist formed from pyroxene dioritic stocks. The greater frequency ratio in NDVI is presented in the range values of 0.2 to 0.3 indicating that the landslides occurred in areas where bare soil or non-healthy low vegetation were present. Moderate values of Topographic Wetness Index (TWI) were present in the area affected by landslides. Areas of coffee crops and secondary vegetation were commonly affected by the landslide history in the study area which was reflected in its FR values.

Performance of MLA

Since the original used dataset presents balancing problems in terms of the variable landslide occurrence (0.34% landslide records, 99.66% non-landslide records), before applying the seven MLA to the training dataset, the SMOTE resampling method was applied to them, creating synthetic samples of the minority class using a k-nearest neighbor algorithm.

On the other hand, according to the results of the Wilcoxon test, the null hypothesis is rejected, concluding the existence of statistically significant differences between all the MLA pairs considered. Once it has been identified, different metrics for each algorithm used to measure their performance were calculated using the test dataset. The validation metrics used, and their results are shown in Table 3. The values of the validation metric highlighted in bold in Table 3, identify the best performing algorithm on that metric. According to these results the best performing algorithm for predicting landslide susceptibility in the study area considering the adopted LCF is the RF algorithm (AUC = 0.95).

Table 3.

Performances Indicators (Metrics) for MLA Considered.

The AUC values for each model are shown in Table 3. This assessment method was used previously by several authors (Hong et al., 2019; Ozer et al., 2020; Rahali, 2019) to check the model performance, and in this work, the ROC curves were also used for the same purpose. The RF method yields the best classification accuracy, and RPART yields the worst performance (AUC = 0.72). All MLA employed in the present study provide acceptable results.

Regarding the results of the other models (Figure 7), the Logit (Figure 7d) and Lasso (Figure 7e) models, which are based on regressions, present similarities in terms of their landslide prediction capacity, showing differences of the order of 10% in all the validation metrics. Regarding the models based on decision trees, XGBoost (Figure 7a), Ranger (Figure 7c), and CTREE (Figure 7f), the spatial distribution of the predicted values of landslide susceptibility do not differ much and show similarities in terms of the spatial patterns of each landslide susceptibility class, except for the RPART model (Figure 7b), in which there is a clear predominance of the “Low” landslide susceptibility class in almost the entire basin that differs spatially with the other models, even in terms of predictive ability using the AUC indicator (Table 3). Compared to the other models, RPART model performed very poorly, indicating that this model is unsuitable for landslide susceptibility mapping in the study area.

Figure 7.

Landslide susceptibility maps using: (a) XGBoost, (b) RPART, (c) Ranger, (d) Logit, (e) Lasso, and (f) CTREE.

According to Huang and Zhao (2018), two main steps should be followed in order to create a landslide susceptibility map: Generate the landslide susceptibility values, and then reclassify them. The method used to reclassify these values depends on the histogram values (Natural breaks, equal intervals, standard deviation, among others). The landslides maps were obtained using the results of the models, which were categorized in five classes: very low, low, moderate, high, and very high classes, using natural break values distribution.

Regarding other studies in the same area with other kinds of methodologies and LCF used, it can be mentioned that Ruiz Vásquez and Aristizábal (2018) obtained an AUC = 0.69 using a multivariate statistical logistic regression analysis. Marin et al. (2021) obtained a performance of the models close to AUC = 0.80 for the total basin, and an AUC = 0.56 for the upper part of the basin using deterministic models (TRIGRS model). Hidalgo and Vega (2021) obtained an AUC = 0.56 using the EPADYM model to calculate the failure probability and factor of safety under seismic and static conditions. Finally, Vega and Hidalgo, (2023) obtained AUC values of 0.95, 0.86 and 0.60 using SVM model, fuzzy gamma model and TRIGRS model respectively, for the landslide-event hazard mapping using the records of May 18, 2015. In other regions with similar topographic conditions with intense rainfall have been carried out studies where the efficiency of MLA was better than the deterministic models, as the case reported in Z. Liu, Gilbert, et al. (2021), where the accuracy was about 82% against accuracies above 90% using RF algorithm.

RF algorithm has better accuracy in prediction compared to linear models. But it cannot be interpreted, so it is often considered as a black-box model. In this work, an interpretable algorithm, SHAP method, is explored for the interpretation of LSA models and the determination of predominant LCF. Figure 8a shows the LSM of the study basin obtained using the RF model. The global prediction results are reasonable to the statistical analysis of different landslide susceptibility classes. The very low susceptibility class covers about more than half of the study area (62.4%; 36.7 km²) and low class 23.2%; (13.6 km²). The areas with moderate and high susceptibility classes cover 7.6% (4.5 km²) and 4.3% (2.5 km²) of the study area respectively. Finally, very high landslide susceptibility class covers 2.4% (1.4 km²) of the basin, almost doubling the area of landslide registered in the basin according to landslide inventory used.

Figure 8.

Landslide susceptibility map using random forest (RF).

Finally, the distribution of landslides in the different susceptibility classes was statistically considered and an analysis of landslide density in the five susceptibility classes was performed. In order to assess whether the generated landslide susceptibility map meets the requirements, two principles are used: First, an ideal landslide susceptibility map has landslide density values that increase from lower to higher susceptibility class (Pradhan & Kim, 2016). Finally, the high-risk areas should account for a small percentage of the total area (Huang & Zhao, 2018). Figure 8b shows the frequency ratio and landslide densities for all classes. Indeed, the very high class has the highest landslide density (LD = 0.08) and landslide frequency ratio (FR = 22.6). Additionally, according to Figure 8c, more than 85% of the events recorded in the inventory occur in this class, according to the values determined with the RF model.

In order to provide global interpretability to the output of RF model, SHAP method was conducted. According to the SHAP values obtained for the test dataset, Elevation and NDVI are the most significant factors for predicting the occurrence of landslides based on SHAP values heatmap (Figure 9). In fact, almost 40% of the dataset has high SHAP values (red color) in both conditioning factors, associated to landslide occurrences while approximately in a 50%, the combination of low SHAP values (blue color) in all LCF is associated to non-occurrence of landslide. The remaining 10%, corresponds to a combination of medium SHAP values of all LCF which lead to landslide occurrence. The obtained SHAP values are consequents with the results of the mean decrease Gini values calculated previously (Figure 10), where both Elevation and NDVI has the higher importance.

Figure 9.

SHAP values heatmap for RF model output.

Figure 10.

Contribution of LCF in RF model output: (a) mean SHAP values and (b) Gini values.

Regarding to the magnitude of each LCF and its influence in RF model output, in Figure 11, it can be noticed that the higher impact on RF model predictions is related to high values of Elevation and low values of NDVI, which is very consistent with the analysis carried out in section 5.2 about spatial distribution of most significant LCF and frequency ratio. High Elevation values in the basin are related to steepest zones prone to landslide occurrence. Low NDVI values are related to a bare soil and poor vegetation areas, very susceptible to infiltration process with the consequent loss of slope stability. The remaining LCF do not present a remarkable impact on RF model predictions.

Figure 11.

SHAP values plot for RF model output.

On the other hand, and as illustration of an example of local interpretability, in Figure 12, it is shown the SHAP values waterfall of a landslide scar pixel (cyan color) in the upper part of study basin, which was correctly classified as landslide with RF model. Again, in this case Elevation and NDVI have the most effect on prediction, whose values are in the conditioning factor class with higher frequency ratio (Table 2). Additionally, the values of Soil Type and Landform in this pixel, correspond to silty clay soils and open steep slopes prone to occurrence of landslides. The remain LCF have no major relevance on RF model output.

Figure 12.

SHAP values waterfall for the pixel N° 8669 classified as landslide for RF model.

Conclusions

Landslides are one of the most widespread and complex natural geodynamics phenomena. In frequently affected regions by landslide, such as tropical mountainous regions, an increase in the number of landslides studies have been shown with the impulse of researchers and regional and local planners. Nowadays, MLA have been applied for landslide susceptibility mapping, making valuable progress. Model algorithm accuracy and diversification have improved susceptibility mapping, particularly with the rapid growth of computer technology and the popularization of GIS techniques.

Relief-F method enabled the detection of potential LCF that could adversely impact the MLA performance, as well as the understanding of LCF that could contribute to the model performance. The LCF representing the terrain characteristics of the study basin were considered very important. Conversely, geologic factors were found to be the least effective factor.

In this study, the landslide susceptibility maps were produced by applying seven different MLA (RF, XGBoost, RPART, Ranger, Logit, Lasso, CTREE). According to some performance metrics, all MLA showed good performance. However, the RF method yields the best classification accuracy with an AUC= 0.95, RPART yields the worst performance with an AUC = 0.72. Based on surface comparison analysis, the most representative results were provided by RF and consequently, RF is the most appropriate approach for landslide susceptibility assessment.

SHAP method provide both global and local interpretability to MLA. According to this method, most influent LCF on landslide susceptibility probabilities were determined. Results revealed that the Elevation and NDVI had the highest positive contribution, while Soil Type and Curvature had no major implications. In addition, it is worth mentioning that terrain and coverage factors had a significant contribution to the landslide occurrences in the study basin, considering that four out of the eight LCF modeled, were the most influential factors on landslide phenomena (Elevation, NDVI, Landcover, and Landforms).

Landslide susceptibility analysis is vital for identifying zones of future landslide occurrence and the proper estimation of landslide-induced risk. In this sense, MLA are efficient and can help to predict disaster risk and decrease disaster costs. The results of a landslide susceptibility analysis using MLA hold immense potential in the field of regional landslide mapping, facilitating the development of effective strategies aimed at minimizing the devastating impacts on human lives, infrastructure, and the natural environment. By leveraging these findings, proactive measures can be devised to safeguard vulnerable areas, mitigate risks, and ensure the safety and well-being of communities. Nevertheless, to fulfill this purpose, interpretation, and explanation of the results of such MLA is required for its correct implementation. Therefore, methods such as SHAP values are used to address this issue in decision making processes based on learning algorithms, providing more abundant and relevant information for landslide risk management.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Sincere thanks to “High Level Training Program for Full-time Professors in their own Doctorates” of the Academic and Research Vice Rector's Offices of the University of Medellín and National Doctorate Program for Teachers of Higher Education Institutions of the Ministry of Science, Technology and Innovation. Data was provided by Research Program “Vulnerability, resilience and risk of communities and supplying basins affected by landslides and avalanches,” code 1118-852-71251, project “Functions for vulnerability assessment due to water shortages by landslides and avalanches: micro-basins of southwest Antioquia,” contract 80740-492-2020 held between Fiduprevisora and the Universidad de Medellín, with resources from the National Financing Fund for science, technology, and innovation, “Francisco José de Caldas.”

© The Author(s) 2023

This article is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License ( https://creativecommons.org/licenses/by-nc/4.0/) which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is attributed as specified on the SAGE and Open Access pages ( https://us.sagepub.com/en-us/nam/open-access-at-sage).