A comparison among fuzzy multi-criteria decision making, bivariate, multivariate and machine learning models in landslide susceptibility mapping

Abstract Landslides are dangerous events which threaten both human life and property. The study aims to analyze the landslide susceptibility (LS) in the Kysuca river basin, Slovakia. For this reason, previous landslide events were analyzed with 16 landslide conditioning factors. Landslide inventory was divided into training (70% of landslide locations) and validating dataset (30% of landslide locations). The heuristic approach of Fuzzy Decision Making Trial and Evaluation Laboratory (FDEMATEL)-Analytic Network Process (ANP) was applied first, followed by bivariate Frequency Ratio (FR), multivariate Logistic Regression (LR), Random Forest Classifier (RFC), Naïve Bayes Classifier (NBC) and Extreme Gradient Boosting (XGBoost), respectively. The results showed that 52.2%, 36.5%, 40.7%, 50.6%, 43.6% and 40.3% of the total basin area had very high to high LS corresponding to FDEMATEL-ANP, FR, LR, RFC, NBC and XGBoost model, respectively. The analysis revealed that RFC was the most accurate model (overall accuracy of 98.3% and AUC of 97.0%). Besides, the heuristic approach of FDEMATEL-ANP model (overall accuracy of 93.8% and AUC of 92.4%) had better prediction capability than bivariate FR (overall accuracy of 86.9% and AUC of 86.1%), multivariate LR (overall accuracy of 90.5% and AUC of 91.2%), machine learning NBC (overall accuracy of 76.3% and AUC of 90.9%) and even deep learning XGBoost (overall accuracy of 92.3% and AUC of 87.1%) models. The study revealed that the FDEMATEL-ANP outweighed the NBC and XGBoost machine learning models, which suggests that heuristic methods should be tested out before directly applying machine learning models.


Introduction
Landslides belong to one of the most significant geodynamic phenomena influencing infrastructure, land use as well as causing loss of lives. Landsliding is conditioned by adverse geological structures and factors causing their activation, such as meteorological factors or inappropriate anthropogenic interventions. Every year, due to climatic factors and inappropriate anthropogenic interventions in Slovakia, landsliding is activated in areas with existing infrastructure or with planned construction (Vlcko et al. 2009).
The total disruption of the territory of Slovakia by landslides is approximately 5% of the total surface area. There is also a growing landslide risk cause of the intensive direction of construction activity, especially in mountainous areas. Most of the landslides are caused by interventions in natural slopes as well as failure to respect the reduced stability of landslide areas (Vojtekov a and Vojtek 2019). In terms of land use, slope deformations contribute almost equally to the violation of agricultural (approximately 50%) and forest land (approximately 47%) and only 3% account for other areas (Ministry of Environment of the Slovak Republic 2020).
Landslides can be divided in terms of the degree of activity into potential (63% out of the territory of Slovakia), stabilized (24.9%), active (11.6%) and combined (0.5%). Moreover, slope deformations within geological formations are related to the Paleogene (60.1%), Neogene (18.7%), Cretaceous (9.9%), Paleozoic (4.2%) and Triassic (2.3%) whilst other geological formations account for 4.8%. In Slovakia, landslides most often occur in Paleogene geological formations. This is also valid for the presented study area -Kysuca River basin, which is mostly formed by flyschoid rocks with a very low degree of permeability (Ministry of Environment of the Slovak Republic 2020).
Landslide susceptibility mapping (LSM) focuses on the identification of the vulnerable areas to potential landslides based on the predispositions which characterize the propensity for landsliding (Fell et al. 2008;Vojtekov a and Vojtek 2020). For that reason, the LS modelling is based on different conditioning factors, which are represented by digital elevation model-derived factors, such as altitude, slope, aspect; meteorological factors, such as annual rainfall or solar radiation; river networkderived factors, such as stream density or distance to rivers; geological factors, such as geological/lithological units and distance to faults; soil factors, land use/land cover factors (Westen et al. 2008). Another important data, either for the training of the models and for the validation of the model results, is the landslide inventory map, which contains the historical landslide locations (Zêzere et al. 2004).
Several groups of methods can be characterized in LSM. The multi-criteria methods (e.g. Analytic Hierarchy Process (AHP), Analytic Network Process (ANP), Weighted Linear Combination (WLC)) which are characterized as knowledge-driven and thus more subjective, were applied in several studies (Neaupane and Piantanakulchai 2006;Akgun et al. 2012;Kayastha et al. 2013;Sezer et al. 2017;Ozer et al. 2020;Yanar et al. 2020). Using the multi-criteria analysis, one landslide conditioning factor directly or indirectly depends on other factors. In particular, DEMATEL method considers the interplay and dependence of decision factors, evaluates the association between the factors, and determines the weight of each relationship. In addition, the ANP method organizes the decision factors into a network of clusters and nodes (Kanani-Sadat et al. 2019;Ali et al. 2021).
On the other hand, the group of more objective methods includes statistical methods (i.e. bivariate and multivariate statistics) (Nohani et al. 2019;Achu et al. 2020), and machine learning methods, which have been emphasized in recent years either as stand-alone or hybrid/ensemble models (Zêzere et al. 2017;Zhang et al. 2018). Recently, Chen and Chen (2021) used bivariate statistical-based kernel logistic regression models with different kernel functions for a case study in China. In another LSM study, Chen et al. (2021aChen et al. ( , 2021bChen et al. ( , 2021c applied the teaching-learning-based optimization and satin bowerbird optimizer algorithms to find optimized parameters of adaptive neuro-fuzzy inference system (ANFIS). In addition, coupled convolutional neural network model with the grey wolf as well as cuckoo optimization algorithms were applied in a study by Li et al. (2021).
In order to construct a landslide probability model, it is assumed that landslide occurrence is dependent on conditioning variables and potential landslides could occur under the similar conditions with historical landslide events. Amongst the bivariate statistical models, frequency ratio (FR) has commonly been used in LS studies (Park et al. 2013;Nourani et al. 2014;Shahabi et al. 2015;Aditian et al. 2018). Furthermore, logistic regression (LR) is popular multivariate statistical model, as evidenced by many LS studies (Choi et al. 2012;Ozdemir and Altural 2013;Youssef 2015). LR is commonly used for classification problems in order to predict solutions based on the concept of probability (Merghadi et al. 2020).
Regarding the machine learning models, the random forest classifier (RFC) is based on a combination of decision tree classifiers and, therefore, it is considered a powerful supervised algorithm for solving binary classification tasks (Breiman 2001). The RFC model has been used in several LS studies, such as Chen et al. (2018aChen et al. ( , 2018b, Sevgen et al. (2019), Nsengiyumva and Valentino (2020), Kocaman et al. (2020) or Zhao et al. (2020). Similarly to RFC model, Naïve Bayes classifier (NBC), a supervised probabilistic algorithm built on Bayes theorem, have been applied in several studies in recent years (Tsangaratos and Ilia 2016;He et al. 2019;Chen et al. 2020aChen et al. , 2020bLee et al. 2020;Lei et al. 2020aLei et al. , 2020b. XGBoost is a scalable tree boosting system built on the basis of gradient boosting in order to reduce computational burden. The ability of XGBoost is to support distributed training and integrate with the cloud dataflow system as well as to utilize more regularized models to resolve the over-fitting problem. The XGBoost model has rarely been applied to LS modelling so far (Sahin 2020;Mirzaei et al. 2021;Pham et al. 2021).
Machine learning methods are considered robust data-driven methods, often with accuracy higher than the multi-criteria methods, bivariate or multivariate statistics . However, this is not the rule and the present study was carried out to investigate this issue as well. In addition, it is found that many of the studies are directly implementing machine learning models without applying any heuristic or statistical methods and concluded that these are the robustness methods, but it would not be universal accepted result because heuristic methods can sometimes outweigh statistical as well as machine learning models.
As for the literature dealing with LS in different regions of Slovakia, several studies applied bivariate and multivariate statistics (Holec et al. 2013; Baran cokov a and Kenderessy 2014; Bu sa et al. 2019), but rarely the machine learning algorithms, such as neural networks (Tornyai et al. 2016). In this regard, a complex comparison of representative methods of the above-mentioned groups was not performed for the study area.
Therefore, comparison of several methods for LSM is the main aim of this study. In particular, the fuzzy multi-criteria decision-making approachesfuzzy DEMATEL-ANP, bivariate statistics -FR, multivariate statistics -LR, machine learning models -RFC, NBC and XGBoost model were applied. The advantage of using several models allowed determining the most suitable one for the study area. The performance and accuracy of the models used were assessed using the Receiver Operating Characteristic (ROC) curves and several other statistical measures.

Study area
Kysuca River basin, particularly its Slovak part, was chosen as the case study covering approximately 989 km 2 . The main river is Kysuca, a length of 66 km, which springs on the northern slopes of the Hri covec peak. The basin is included in the NUTS II Central Slovakia and NUTS III Zilina Region. A total of 55 municipalities extend to the Kysuca River basin.
The geological structure of the study area consists of the following rocks: claystone-limestone rocks (0.7%), conglomerate-sandstone rocks (7.3%), deluvial deposits (13.6%), floodplain deposits (3.6%), flyschoid rocks (73.2%), limestone-dolomite rocks (1.5%) and river terrace deposits (0.1%). The geomorphological units, which extend to the basin, can be divided to mountains, foothills and basin (Maz ur and Lukni s 1986). Rycierova hora (1226 m), which is situated in the eastern part of the basin, represents the highest peak. On the contrary, the lowest elevation (approximately 326 m) is situated in the place where Kysuca River flows into the V ah River.
The study area is located in the temperate climatic zone. Based on the climatic classification of Slovakia, as presented by Lapin et al. (2002), the basin is included in two climatic regions. Lower lying areas belong to a moderately warm region with the July mean temperature of 16 C or more and an average annual number of summer days of less than 50 days. The rest of the basin is included in a cool region with the July mean temperature of less than 16 C. The average annual rainfall rises from the basin elevations with 800 mm of rainfall to the mountainous elevations with 1500 mm of rainfall (Bochn ı cek et al. 2015).
Regarding the occurrence of landslides, the Kysuca basin belongs to highly vulnerable region. It is documented by many landslides, which occurred in the study area, such as in 2019 in the town of Cadca, in 2013 in the municipality of Rudinka, in 2010 in the municipality of Svr cinovec, in 2009 in the municipalities of Lodno, Klubina and the town of Cadca and so on. The disruption of the Zilina Region, where the study area belongs, by landslides is the highest (12%), compared to other regions of Slovakia. The study area is included in three districts, where the disruption by landslides is 10%, 20% and 22% in case of the Zilina District, Cadca District and Kysuck e Nov e Mesto District, respectively (Ministry of Environment of the Slovak Republic 2020).

Landslide inventory map
The past landslides data is critical for LSM, especially, for the training as well as testing phase. The historical landslide locations (i.e. inventory map) were taken from the Atlas of Slope Stability Maps at a scale of 1:50,000 ( Simekov a et al. 2006). All landslides were recorded in a polygon format. The study area contains 180 polygon landslides, which differ in size from 142 m 2 to 497,840 m 2 . In this study, 70% of landslide locations were utilized in the training and the remaining 30% were reserved for validation (Chen et al. 2018a(Chen et al. , 2018bAli et al. 2020). A recent study showed that model performance using training/validating dataset with a combination of 70%/ 30% offers best result in comparison to combinations of 60%/40, 80%/20% and 90%/ 10% (Shirzadi et al. 2018). It should be noted that the landslides are continuously incorporated into the atlas and it also contains landslides which occurred after the year of its publication. Landslides have been identified based on the fieldwork by various organizations, such as private companies and universities.
In study area, the landslide inventory data consisted of slow-moving, which were caused by the change in the hydrologic-geologic processes in the area of the landslide, especially, an increase in levels of groundwater. The share of deep-seated landslides out of all landslide locations was 16.1%. On the other hand, shallow-rapid landslides, which were caused predominantly by intense rainfall, but also lateral/vertical erosion or their combination, accounted for 83.9% of all landslide locations. In this cases, sudden saturation loosened the soil and triggered the landslide, i.e. gravity overtook the hillside and the muddy soil mass broke loose. In addition, 93.3% of landslides occurred on Paleogene rocks, 3.4% on Cretaceous rocks, and 3.3% on Jurassic rocks. The geological structure of 67.8% of landslides was represented by mixed and talus soils/eluvium whilst the geology of 32.2% of landslides was created by clayey and fragmentary unconsolidated rocks (consistent and frictional soils).
For the purposes of LS modelling using the raster analysis, the areal landslides were converted to 1000 landslide points using ArcGIS software. To improve the training and testing phase of LS modelling, another layer of 1000 non-landslide points was also created. The non-landslide locations were identified only for areas with a slope degree lower than 3 . As a result, the LS analysis was performed with the use of 2000 points ( Figure 2).

Landslide conditioning factors
The landslide conditioning factors represent crucial input data for LSM. Selection of relevant conditioning factors is usually carried out in a subjective manner and there are no general principles for the selection of these factors (Vojtekov a and Vojtek 2020). The process of selecting conditioning factors may also depend on the availability of quality and accurate source data as well as on the nature of landslide activities (Vakhshoori et al. 2019).
As for the current study, 16 conditioning factors were considered to be relevant for LSM in this study area, which sources and primary data are shown in Table 1.
Morphometric/hydrographic factors are represented by altitude, slope angle, slope aspect, general curvature, plan curvature, profile curvature, sediment transport index (STI), topographic wetness index (TWI), stream power index (SPI) and distance to rivers ( Figure 3). Climatic factors are the average annual precipitation and annual solar radiation (Figure 4). Geological factors include lithological units and distance to faults ( Figure 5). The last group contains land cover and distance to roads ( Figure 6). All conditioning factors were processed as raster with 20 Â 20 m spatial resolution using ArcGIS software.
The basis for creating factors of elevation, slope angle, slope aspect, general/profile/plan curvature and annual solar radiation was the digital elevation model (DEM). By using the topographic map with a scale of 1:50,000, the DEM was interpolated from elevation points and contour lines. STI, TWI and SPI were calculated in ArcGIS software using Eqs. (1) where A s is the slope contributing area and b is the slope gradient. Distance to rivers and roads was computed by using ArcGIS software based on the river and road network, respectively, in the topographic map with a scale of 1:50,000. The rainfall factor, i.e. average annual total precipitation, was prepared based on the Climate Atlas of Slovakia (Bochn ı cek et al. 2015). The mean annual rainfall was used as a conditioning factor since it represents the climatological conditions of the study area. In turn, it cannot be considered a triggering factor expressed, for example, by sub-hourly, hourly or daily precipitation or rainfall intensity and duration. The use of annual total rainfall could provide information about the potential of landsliding in different climatological conditions (Tsangaratos and Ilia 2016). Lithological units as well as the distribution of faults were taken from the Engineering-geological Zoning Map

Methods
In this study, FDEMATEL-ANP, RFC, NBC, XGBoost, FR and LR models were used for comparative assessment of LS. First, the following tasks were performed: preparing the landslide inventory map and determining landslide conditioning factors, normalizing the factors, extracting weights using the FDEMATEL-ANP approach and preparing the LS models. Second, feature selection to choose effective conditioning factors for model building, then generating LSM using the aforementioned models. The detailed of methodology adopted in this study is shown by the following graph ( Figure 7).

Multicollinearity test
The multicollinearity problem between more than two parameters leads to the difficulties in predicting landslide occurrence, usually, causing errors in results. Therefore, suitable selection of these parameters should be done to confirm that they are independent of each other (Achour et al. 2018). As indicators of multicollinearity, the tolerance (T) and variance inflation (VIF) indices were used to perform this test in this study. A VIF > 5 shows high multicollinearity between conditioning factors. In addition, T > 0.2 exhibits low multicollinearity (Pham et al. 2016). The results of this test are summarized in Table 2, which displays that four factors, namely general curvature, STI, annual solar radiation and TWI, were rejected due to high VIF and low T value. On the other hand, other 12 factors, including lithology, elevation, slope angle, slope aspect, plan curvature, profile curvature, distance to faults, rainfall, stream power index (SPI), distance to roads, distance to river and land cover, were in the safe range of collinearity and thus they were used for LS modelling.

Fuzzy decision-making trial and evaluation laboratory (FDEMATEL)
The DEMATEL belongs to one of the multi-criteria decision making (MCDM) techniques which are used to determine the inter-relationship and dependence of selected criteria. Using this method, the direct relation matrix is constructed to find the dependencies and intensity of direct relation amongst the criteria at a comparison scale from 0 to 4 (Yang et al. 2008). For example, if the criterion 2 is ranked with 1 and criterion 7 with 4, it means that the criterion 7 has a very high influence upon the criteion 2 . However, to remove the subjectivity in expert judgments and inconsistency in the variables, fuzzy logic is used to create the fuzzy direct relation matrix (Zadeh 1965). Table 3 presents the linguistic variables and triangular fuzzy numbers, which are needed for calculating the fuzzy direct relation matrix. The FDEMATEL method can be briefly represented in the following main steps. A detailed description of this method can be found, for example, in the studies by Ali et al. (2021) or Kanani-Sadat et al. (2019).
Step 1: Constructing the fuzzy direct relation matrix using several pairwise comparisons performed by experts.
Step 2: Calculating the normalized fuzzy direct relation matrix.
Step 3: Computing the fuzzy total-relation matrix which shows the direct and indirect relations amongst the factors.
Step 4: Calculating the fuzzy values for the sum of columns and rows.
Step 5: Determining the causal relationship amongst the criteria for analyzing the comprehensive inter-relationship of decision criteria and to judge which are the most influential conditioning factors and how they affect other factors.
Upon preparing the fuzzy total-relation matrix, the direct and indirect relations were found amongst the factors that were determined by calculating the summation of columns and rows, symbolically known asD f andR f : Thus,R f þD f andR f ÀD f were estimated to understand the causal relation and influence of a factor with respect to others. The higher and lower value ofR f þD f implies the greater and lower priority respectively of a conditioning factor in landslide incidences. Supermatrix, the first step of the analytic network process (ANP) was determined based on the priority of factors derived from the last stage of fuzzy DEMATEL to make ensemble FDEMATEL-ANP.

Analytic network process (ANP)
The ANP is another multi-criteria technique, which divides the criteria into a network of clusters and nodes due to the presupposition that there are inter-dependencies amongst the criteria (Saaty 1996). For this study, the ANP technique was applied to calculate the final weights. The following three main steps were undertaken when applying the ANP. The ANP technique is further described, for example, in the works of Ali et al. (2020) or Kanani-Sadat et al. (2019).
Step 1: Preparing the super-matrix with the use of pairwise comparisons of the criteria and applying a continuous scale which varies from 1 to 9 (from equal importance to extreme importance) (Saaty 1996).
Step 2: Calculating the weighted super-matrix based on the summation value of each row of the threshold value.
Step 3: Converging the weighted super-matrix to a stable matrix by its limitation powered to a large enough number.
After calculating the final weights of selected factors by stabling the matrix as stated in step 3, the final LS map was generated based on the derived weights.

Frequency ratio (FR) model
FR is applied to estimate the probability of landslide occurrence (Choi et al. 2012). The FR presents the observed spatial relationships between the landslide occurrences and causal landslide factors (Lee and Min 2001;Lee and Talib 2005). For each class of factor, the frequency ratio for was calculated using Eq. (4) (Ali et al. 2019). Then, the landslide susceptibility index (LSI) was calculated by summation of frequency ratios using Eq. (5): where s j i j is the number of landslide pixels in class j of factor i; P j s j i is the total number of landslide pixels; N j i is the number of pixels of class j in factor i; P j N j i is the total number of pixels of the whole area.
The landslide occurrences and causal landslide factors have a higher correlation when the value of Fr is higher than 1. The landslide occurrences and causal landslide factors has a lower correlation between landslide occurrences and their causal factors when the value is lower than 1.

Logistic regression (LR) model
LR, i.e. a multivariate regression analysis, is applied for the presence or absence prediction of a characteristic or phenomena based on values of a set of predictor variables (Achu et al. 2020). The LR model presents the spatial regression relationship between landslide-presence/absences and factors influencing landslides. In this study, the dependent variable (landslides) is a binary variable representing the presence or absence of landslides and the values must be input as either 0 or 1. Logistic regression coefficients are used to estimate the ratio between dependent and independent variables in the model. The relation of the probability that is possibility of landslide occurrence and the independent variables (contributing factors) is presented in Eq. (6): where p is the probability of the landslide event occurring and its value ranges between 0 and 1. z is the linear combination and calculated by Eq. (7): are the slope coefficients of the LR model.

Random forest classifier (RFC)
RFC is applied for the multivariate classification of large datasets using an ensemble of decision trees, thus ensuring the model stability (Breiman 2001;Catani et al. 2013). RFC uses a combination of a large number of decision trees, grown on a bootstrap sample, to compute a classification by resampling the dataset with substitution and randomly altering the predictor sets over different tree induction developments (Achour and Pourghasemi 2020). The RFC requires the number of trees and variables used to divide the nodes whilst each node is divided by the best of a random subset of predictors. The so-called out-of-bag (OOB) error rate is computed for each tree with the use of observations left out of the bootstrap sample. Finally, based on the majority vote of all trees, the class membership and output are determined (Park and Kim 2019). The RFC selects a higher number of nodes to its output class. Therefore, for an input data 'i' is computed its output 'o' from maximum ensemble which is shown in Eq. (8): where f ðdÞ is an indicator function defined as: where 0 YES 0 and 0 NO 0 are the landslide and non-landslide location, respectively.

Naï ve bayes classifier (NBC)
NBC is a machine learning method based on the Bayes' theorem, which assumes the condition that all the training variables are independent and each variable contributes equally to the problem of classification (Tsangaratos and Ilia 2016;Lee et al. 2020).
The NBC use the simple iterative parameter estimation schemes, therefore, this method develops easily. On the other hand, NBC is robust to noise and irrelevant attributes (Tien Bui et al. 2012). According to Khosravi et al. (2019), the main steps in applying the NBC are (i) selection of training variables, (ii) calculation of the posterior probability of each class (landslide/non-landslide locations), (iii) estimation of class level as well as a covariance matrix and (iv) creation of a discriminant function for each class level. The classification probability is thus obtained for all classes. The classification rule in naïve Bayes classifier is represented with respect to 'n' element which is shown in Eq. (9): where p i is the landslide causative factor (i ¼ 1, 2, … , n), q j is the output class (j¼landslide/non-landslide). Naïve Bayes evaluates the landslide susceptibility L p i q j ð Þ for output classes. The evaluation is based on the class with the largest posterior susceptibility whilst the prior susceptibility Lðp i Þ can be estimated by the fraction of the observations with output class q j in the training dataset.

Extreme gradient boosting model
Based on the theory of Gradient Boosting Machine (GBM) (Friedman 2002), the XGBoost algorithm was improved and presented as a powerful decision tree by Chen et al. (2015). It integrates several decision tree models to generate a classification, with higher accuracy, based on iteration process. XGBoost has reached expected results in LSM (Zhang et al., 2019). The function of XGBoost algorithm is given in Eq. (10) (Zhang et al., 2019): where i signifies the ith sample, ŷ i (tÀ1) serves as the expected value of the (t À 1)th model for the sample i, f t (xi) means the newly added tth model, X(f t ) refers to the regular term, C represents some constant terms, and the outermost L() is the error.
In this study, The R V R software version 3.6.2 was used to apply XGBoost. The following XGBoost parameters were checked and ameliorated: number of iterations (m), learning rate (g), maximum tree depth (h), minimum loss reduction (c), minimum sum of instance weight (f), subsample percentage (l) and subsample ratio of columns (⑁) ( Table 4).

Model performance
The ROC curves were applied to evaluate the accuracy of the LS models. Using the area under the curve (AUC), which is the main statistical indicator of the ROC curve.
In the present study, sensitivity is the ratio of true positive points to all landslide classes (Eq. (11)), whereas 1 À specificity defines the true negative points to all nonlandslide classes (Eq. (12)) (Gajowniczek et al., 2014): where TP and TN are the true positive and true negative, respectively, which represents number of landslide pixels correctly classified. FN and FP are the false negative and false positive, and respectively, representing number of landslide pixels incorrectly classified. Moreover, if landslide is represented as 1 and non-landslide is represented as 0, the AUC can be expressed as Eq. (13): where t is the threshold value.
The model is considered inaccurate if the ROC curve is lower than 0.5 (Ali et al. 2020;Chen et al. 2020aChen et al. , 2020b. The model is more accurate performance if the ROC curve is closer to 1 . Along with ROC curve, the model performance was also tested with overall accuracy (OC). The overall accuracy is an important statistical measure used for model performance and validation (Ali et al. 2020). The mathematical inference of OC is expressed in Eq. (14) which considers both sensitivity (true positive rate) and specificity (false positive rate).
Regarding the validation of the LS model via FR model, the success rate curve was used. It is an estimation of the success of a model that indicates how well the model matches the prior events (Lei et al. 2020a(Lei et al. , 2020bChen et al. 2021aChen et al. , 2021bChen et al. , 2021c. Success rate curve was created by using training data and LS map (Chen et al. 2021a(Chen et al. , 2021b(Chen et al. , 2021c. To produce this curve, the computed index values of all cells were organized in descending order and were split into 256 equal classes going from high susceptibility zones to non-susceptible one. Thereafter, we plot the susceptible classes starting from the highest to the lowest values on the X-axis and the cumulative percentage of landslides occurrence on the Y-axis.

Results and analysis
The LS modelling was studied for the Kysuca basin, Slovakia using GIS-based fuzzy multi-criteria decision-making approach, machine learning algorithms, bivariate statistics and logistic function. Initially, 16 relevant affecting factors of the landslide were selected. However, the multicollinearity diagnostics of these factors lead to the elimination of four unnecessary factors (general curvature, sediment transport index, annual solar radiation and topographic wetness index) for effective LS modelling (LSM) ( Table 2). The multicollinearity test showed that profile curvature (T ¼ 0.98), distance to river (T ¼ 0.96), elevation (T ¼ 0.86) and rainfall (T ¼ 0.83) were the most significant factors for LS modelling having higher T value. All other conditioning factors have also significant role for LS modelling with T value >0.20, except for slope (T ¼ 0.46), SPI (T ¼ 0.46) and distance to faults (T ¼ 0.39). As a result, the spatial distribution of LS was estimated using these 12 conditioning factors and six presented models, followed by comparative analysis and accuracy assessment. The final results of each model are presented in the sub-sections below.

Landslide susceptibility mapping via fuzzy DEMATEL-ANP
The fuzzy DEMATEL-ANP calculates the inter-dependences and inter-relationship of each criterion in comparison to others. The R-value and D of the defuzzified totalrelation matrix identified the role of all affecting factors in landslide susceptible map ( Table 5). The higher defuzzified value of R df þD df indicated the higher relation with other affecting factors whereas the nature of relationship and dependencies were assessed by the values of R df À D df . If the value of R df À D df is in plus (þ), it comes to cause group which affects the other factors. Concomitantly, if the value of R df À D df is in minus (À), it will consider to be effected group that influenced by other factors. In this analysis, out of total landslide conditioning factors, the aspect (0.031), distance to Based on the FDEMATEL-ANP, the normalized weights, as shown in Table 6, were multiplied with respective conditioning factors to produce the LS map. The final susceptibility map, prepared using the FDEMATEL-ANP (Figure 8(a)) model, was categorized into 5 landslide susceptible classes using the natural break method: very high, high, moderate, low and very low. Figure 8(a) shows that 23.41% out of the whole study area was located in very high landslide susceptible areas. Furthermore, 28.83% and 22.73% were located in high and moderate susceptible areas, respectively, whereas 16.66% and 8.34% were located in low and very low landslide susceptible areas, respectively.

Landslide susceptibility mapping via FR
The FR values were computed using the MS Excel spreadsheet. The FR model was used to identify the correlation between landslide points and each affecting factor. Table 7 shows the ratio value of each selected landslide conditioning factors. For FR, a value close to 0 or perfectly 0 indicates low to very low LS whilst its value higher than 1 indicates a higher susceptibility to the landslide. Therefore, a greater value of FR signifies higher susceptibility to landslide and vice versa.
The relationship between landslide occurrence and lithology shows that flyschoid rock and limestone-dolomite rock are the major landslide occurring lithological group with FR value 0.96 and 0.93, respectively. The FR values of elevation were found greater for elevation >1100 m and 900-1100 m with FR values of 2.55 and 2.16, respectively. It was thus evidenced that the landslide occurrence increases with increasing elevation. For the Kysuca river basin, the elevation ranges from 326.60 m to 1238.12 m and the FR values for classes <500 m to >1100 m was calculated greater than or close to 1, indicating high LS. The relationship between slope angle and landslide incidence was also shown using FR. The slope angle classes of > 30 and 20-30 have the highest values of FR with weights of 3.38 and 1.14, respectively. The result also indicates, as in case of elevation, that the LS increases with an increasing slope angle. For the slope aspect, the slope facing east (FR ¼ 1.17), southeast (FR ¼ 1.10), west (FR ¼ 1.09) and southwest (FR ¼ 1.03) had the highest FR values. Hence, these aspect classes are highly susceptible to landslide whereas the flat and northwest facing classes had the lowest values with FR ¼ 0 and FR ¼ 0.88, respectively, indicating lower susceptibility to the landslide.  Curvature is one of important parameters for LS modelling. Both the convex and concave curvatures are more susceptible to landslide than flat curvature. The present analysis shows that for plan curvature, the FR values of convex and concave curvature were, respectively, 1.32 and 0.93. For profile curvature, the FR values of concave and convex curvature were 1.91 and 1.11, respectively, whereas the FR value of flat curvature was <1 (Table 7). For distance to faults, the classes of <1000, 1000-2000 and 2000-3000 m had greater FR weights with the values of 1.29, 1.23 and 1.21, respectively, implying very high LS. Consequently, higher distance to faults, i.e. 3000-4000, 4000-5000 and >5000 m with <1 FR values, denotes a lower susceptibility to the landslide. The results of frequency ratio and rainfall data demonstrated a positive correlation between them which implies that the increasing landslide incidences were directly related to increasing rainfall. The FR values of the rainfall class 900 À 1000 mm (FR ¼ 1.26) was >1 whereas the lowest FR value (0) was found for <800 mm rainfall class.
The relationship between SPI and landslide occurrence demonstrated that classes of < À2, and À2 À 0 had high FR weights with values of 1.40 and 1.87 whilst increasing the value of SPI denoted lower FR values ( Table 7). The likelihood of landslide occurrence increases with decreasing distance from roads and rivers. Therefore, areas close to roads and rivers are highly prone to landslide. The results in the present study demonstrated that in case of distance to roads, the FR weights 1.29 and 1.11 were found in <50 m and 50-100 m classes, respectively, whereas lower FR weights 0.77 was in >600 m class. Same as for the distance to rivers, the frequency ratio values 1.26 and 1.01 were found in <100 m and 100-200 m classes individually and lower value 0.51 was found for the class >500 m from the river. This means the occurrence of landslides decreases with increasing distance from roads and rivers. In case of land cover, higher weights of FR were noticed for commercial units, transitional woodland-shrub and discontinuous urban fabric with FR weights of 1.73, 1.34 and 1.27, respectively. Subsequently, land cover classes like sport and leisure facilities and permanent croplands had FR weights of 0, which indicates the least probability of landslide occurrence (Table 7).
Based on FR, the LS map was produced and the spatial distribution is shown in Figure 8(b). The final susceptibility map prepared using the FR model was classified into five landslide susceptible classes using the natural break method: very high, high, moderate, low and very low. Figure 8(b) shows that out of the total area, 11.83% were located in very high landslide susceptible areas. Moreover, 24.63% and 27.87% were located in high and moderate susceptible areas, respectively, whereas 23.50% and 12.16% were located in low and very low landslide susceptible areas, respectively.

Landslide susceptibility mapping via LR
In the LR model, the beta (b) coefficients of each conditioning factor were used for calculating the LSI. The formula can be rewritten in Eq. (15) The aforementioned b coefficients were multiplied by the respective conditioning factors using the raster calculator in ArcGIS software to obtain the final landslide susceptibility map (Figure 8(c)). From Eq. (1), the results of LR coefficients showed that slope angle, slope aspect, plan curvature, SPI, rainfall, distance to faults, distance to roads, distance to rivers, and elevation had positive b value indicating that these factors have a prominent role in the LS modelling in the Kysuca basin. It is also seen that slope angle had the highest b coefficient (0.856), followed by SPI (0.478), distance to roads (0.245), slope aspect (0.222), plan curvature (0.208), elevation (0.112), distance to rivers (0.09), rainfall (0.08) and distance to faults (0.014), respectively. Land cover and profile curvature had a negative b coefficient and, therefore, these factors had a negative effect and were less important for LS.
Similar to the other LSI models, the final susceptibility map was divided into five classes using the natural break method: very high, high, moderate, low and very low. The classified results showed that out of the total area of the Kysuca basin, 14.30% and 26.41% belonged to a very high and high susceptible class, respectively, whilst 27.36%, 21.07% and 10.84% belonged to moderate, low, and very low susceptible classes, respectively (Figure 8(c)).

Landslide susceptibility mapping via RFC
The random forest model was prepared using 700 landslide pixels as the training data. The 'k' selected after hyperparameters tuning was 1000, and the number of features tested at every node was chosen to be 3, which means that three factors were tested at each spilt. The aggregate OOB, which estimates the prediction error of random forest, is shown in Figure 9. The lower value of OOB indicates higher accuracy of models. For the present study, the OOB was <2% indicating that the accuracy of RFC was >98% which is a very good and acceptable model. The importance of landslide conditioning factors assessing through RFC is shown in Figure 10. The results of this figure revealed that slope angle, distance to roads, plan curvature and elevation are the essential factors for landslide susceptibility according to both mean decrease accuracy and mean decrease Gini. The other factors like slope aspect, lithology, profile curvature, distance to faults, distance to rivers, SPI, rainfall and land cover are moderate to very low significant factors, respectively.
After completing the calculation, the relevant data obtained from the RFC model was transferred in the Geoprocessing tool of ArcGIS software to prepare the final LSI. The final LS map is shown in Figure 8(d). It was further classified using natural breaks classification procedure into very high, high, moderate, low and very low susceptible classes. The results of LS modelling via RFC revealed that 20.34% and 30.26% of the case study had very high and high LS. Moderate, low and very low susceptible areas were found to have 24.12%, 17.05% and 8.21%, respectively.

Landslide susceptibility mapping via NBC
For the NBC model, the weight of attribute rank was used for preparing the LS map, which is shown in Table 6. In case of NBC, the slope angle was found as the important conditioning factor for LS with a weight of 0.9636. The other factors like SPI (0.8471), elevation (0.8356) and distance to faults (0.7724) were also significant in LS with the higher weight of ranked attribute. On the other hand, distance to roads, slope aspect, distance to rivers were the least important factors in LS in this case study with weights of 0.0350, 0.0758 and 0.0955, respectively (Table 6).
Using the NBC, strong assumptions were taken between landslide training points and feature classes of conditioning factors. It was predicted that probabilities of training datasets, i.e. landslide points, belong to each class of landslide conditioning factor and the classes with the highest probabilities are considered the most probable classes to the landslide.
The weights resulting through NBC were multiplied with respective conditioning factors to construct the LS map (Figure 8(e)). The final susceptibility map, prepared using the NBC model, was reclassified again into five susceptible classes: very high, high, moderate, low and very low. The results revealed that out of the total area of the Kysuca river basin, 15.68% was located in the very high LS class. Further, 27.89% and 26.52% were belong to high and moderate susceptibility classes, respectively, whilst 19.43% and 10.45% were located in low and very low LS classes, respectively.

Landslide susceptibility mapping via XGBoost
XGBoost is an ensemble machine learning method that is used to improve the accuracy of models and reduce over-fitting problems. XGBoost is used both for classification and regression. In the present study, the training datasets for running the XGBoost model were prepared using the twelve landslide conditioning factors to shuffle data and separate factors to reach our objective and predict LS (0 ¼ non-landslide; 1 ¼ landslide). The XGBoost model sums its all trees for predicting the weights, in which the number of trees is controlled by the number of predictors. The number Figure 10. The importance of conditioning factors as determined by the RFC model. of trees was 100 by defaults for model building in the present study. Each tree in XGBoost classification may not produce a good prediction, but the aggregation of all trees produces a robust prediction. Therefore, the classification accuracy in XGBoost was carried out by checking with 10-fold cross-validation. Hyperparameters tuning was applied to optimizing the values of the parameters. The hyperparameters tuning will not be effective if it is closed to its maximum performance level with the datasets. The ideal XGBoost parameter values after hyperparameters tuning are presented in Table 4. The results proved that the model with its maximum performance comes to have 92.3% accuracy.
To estimate the importance of the factors, the weight was used and LSI was calculated. The final LS map using XGBoost is shown in Figure 8(f). For spatial distribution, the final susceptibility map was classified into five categories. The spatial results of LS via XGBoost showed that 40.31% of the total area had very high to high LS. Moreover, 19.74% had moderate susceptibility whereas low and very low susceptible areas were found to have 20.79% and 19.13%, respectively.
The present study analyzed and compared six heterogeneous models for LSI. For each model, the training and testing datasets were constructed using the common sampling strategy (Ali et al., 2021;Fang et al., 2021). Note that the estimated weight of each landslide model is significant for final LSI. In this study, different weighting methods and cross-validation procedures were used. Particularly, the R-value and D of the defuzzified total-relation matrix were used for fuzzy DEMATEL-ANP which identifies the role of each conditioning factor in LSI. The ratio between landslide pixels and affecting factor was used in frequency ratio (FR) model to determine the spatial relationship between occurrence of landslides and conditioning factors which improves the classification accuracy. In case of machine learning models, such as in LR, the beta (b) coefficients of each factor were taken to calculate the LSI which helps to deduce the over-fitting problem whilst constructing the model; in random forest classifier, the hyperparameters tuning was performed to select 'k' value and determine importance of conditioning factors; in case of NBC, the attribute rank was derived to weight the conditioning factors; whereas hyperparameters tuning was used in XGBoost to optimize the selected factors after 10-fold cross-validation.
From the critical point of view, it has been found in the above discussion that FR is an effective model to apply, but it is a classical approach used in LS. Fuzzy DEMATEL-ANP offers good accuracy, but it is a complex MCDM approach which involves experts that need more time to process. LR is a simple machine learning model to apply, but it needs independent variables, whereas RFC, NBC and XGBoost models are suitable for multi-dimensional data analysis without over-fitting problem, but these are more complex in term of model construction and application.
So, this study aimed to apply the aforementioned heterogeneous models in order to determine the more efficient and suitable models for LSM. The landslide indices were classified into five susceptible zones as shown in Figure 8. Although, all applied models are heterogeneous by nature, but similar spatial results have been found from each model, i.e. very high landslide susceptible areas are located in the north-eastern part, whereas very low susceptible areas are located in the northern and southern parts of the site with lower elevation (<700 m) and lower slope (<10 ).

Comparison of model performance and validation
Six different models, namely FDEMATEL-ANP, FR, LR, RFC, NBC and XGBoost were constructed for spatial prediction of LS. The models were validated with state-of-theart statistical metrics (e.g. specificity, sensitivity and overall accuracy) whilst ROC-AUC curves were computed. Table 8 shows the results of the specificity, sensitivity and accuracy of all models. The results revealed that the highest specificity was fitted to the RFC model (97.5%), which indicated that 97.5% of the non-landslide pixels were correctly classified. The next model concerning higher specificity was FDEMATEL-ANP (94.8%) followed by NBC (92.7%), XGBoost (91.2%), LR (88.6%) and FR (87.3%), respectively. Based on sensitivity measures, the RFC model had the highest sensitivity (99.1%), implying that 99.1% of the landslide pixels were matched with the landslide class, followed by the FDEMATEL-ANP (97.2%), NBC (94.6%), XGBoost (93.4%), LR (92.4%) and FR (90.3), respectively. The overall accuracy (OA) was also shown in Table  8, which illustrated that the RFC model with value 98.3% had the highest accuracy, followed by the FDEMATEL-ANP (93.8%), XGBoost (92.3%), LR (90.5%), FR (86.9) and NBC (76.3%), respectively. The results produced from each of six models cover LS pixels differently, which is also shown for visual comparison in Figure 11.
The results of all six models applied in LS indicated that the models offered a very good and accepted level of performance for both spatial prediction and distribution of landslides in the Kysuca basin. Though, the RFC model had the highest performance rate with respect to specificity, sensitivity and overall accuracy.
As for the ROC curves, they were drawn by comparing the training datasets and final landslide susceptible maps, which is shown in Figure 12. As per the AUC plot, the RFC model with AUC 97% had the highest accuracy outperforming the other models. The next highly accurate model was FDEMATEL-ANP (AUC ¼ 92.4%), followed by LR (AUC ¼ 91.2%), NBC (AUC ¼ 90.9%), XGBoost (AUC ¼ 87.1%) and FR (AUC ¼ 86.1%). The results proved that the RFC model had the highest accuracy in LSM in the Kysuca basin, but other models of FDEMATEL-ANP, NBC, LR, XGBoost and FR had also good accuracy in overall prediction of LS. It was also found from the present analysis that the fuzzy DEMATEL-ANP model had greater accuracy than bivariate (FR) and multivariate (LR) approaches and even deep learning (XGboost) model. Therefore, it is highly recommended to implement selected heuristic methods before moving to multivariate and machine learning models.

Discussion
The present study evaluated different GIS-based models, particularly FDEMATEL-ANP, FR, LR, RFC, NBC and XGBoost for LSM in the Kysuca basin. First, 16 evaluation factors were selected, but after the multi-collinearity test, 12 factors were selected for analysis and others were eliminated because the values of Tolerance and VIF were out of range. The twelve conditioning landslide factors including lithology, elevation, slope gradient, slope aspect, plan curvature, profile curvature, distance to faults, rainfall, SPI, distance to roads, distance to rivers and land cover were considered for the LS assessment. Factors of elevation, slope angle, slope aspect, plan curvature, SPI and distance to faults had a more significant role in LS that was evidenced from all six models applied in the present study. For the RFC model, the factors of slope angle, distance to roads, plan curvature and elevation with the higher value of mean decrease accuracy and mean decrease Gini were positively correlated with the occurrence of the landslide. For the LR model, the result revealed that the slope angle with a weight value of 0.856 was the most significant factor which had a positive relationship with the landslide incidence in the study area. The other factors like SPI (weight value ¼ 0.478), distance to roads (weight value ¼ 0.245), and plan curvature (weight value ¼ 0.208) were also important after slope gradient. The findings based on FR model showed that the slope angle with ratio value of 3.38 was the most important factor, followed by lithology (ratio value ¼ 3.07), rainfall (ratio value ¼ 2.77), elevation (ratio value ¼ 2.55), SPI (ratio value ¼ 1.87) and distance to roads (ratio value ¼ 1.29) which had also positive correlation with landslide occurrence. According to the results obtained from FDEMATEL-ANP model, the slope angle was a more effective factor for LS with the normalized weight of 0.1707 than other factors like distance to faults, distance to rivers, elevation, plan curvature and rainfall with the normalized weight of 0.0877, 0.0679, 0.0524, 0.0436 and 0.0436, respectively. In case of the NBC model, the highest NBC weight values were found for slope angle (0.9636), followed by SPI (0.8471), elevation (0.8356), distance to faults (0.7714), rainfall (0.4860), land From the above discussion, it was found that the slope angle had the highest weight value in all of the six models. Hence, it was the most influential factor for LS prediction in the study area, which has an important impact on the equilibrium of slopes and it is responsible for its stability and instability. The probabilities of landslides occurrence are increased when slopes are steep and unstable (Achour and Pourghasemi 2020). Along with the slope angle in all models, other topographic factors like elevation, SPI and plan curvature had also significant role. As a result, all applied models evidenced that topographic reasons, like a higher degree of slope deformation, high elevation, positive values of plan curvature and SPI (i.e. the relationship between landslide and flow erosion), were the most essential factors for LS in the basin.
Recently, several machine learning algorithms and modelling techniques were used for preparing and assessing LS (Chen et al. 2018a(Chen et al. , 2018bKavzoglu et al. 2018). The ensemble learning methods in LS were prioritized because of their high prediction accuracy (Hong et al. 2018;Vakhshoori et al. 2019;Dang et al. 2020;Tien Bui et al. 2020). A recent study compared RFC, SVM and XGBoost for multi-geohazards susceptibility in Jiuzhaigou, China and found that XGBoost had the best prediction capability with highest AUC of 93% (Cao et al. 2020). Similarly, three different machine learning techniques, RF, SVM and boosted regression tree (BRT) were applied for assessing the accuracy of LS along A1 Highway of Algeria and found that RF had higher predictive accuracy than SVM and BRT models with AUC of 97.2% (Achour and Pourghasemi 2020). Another study by Park and Kim (2019) compared between the RFmodel and the BRT model for LS at Woomyeon Mountain of South Korea and found that RF model had better accuracy in comparison to BRT. Similarly, Chen et al. (2017) compared three models, namely, logistic model tree, RF and classification and regression tree models. Their finding revealed that RF model produced the highest accuracy in comparison to the other two models. Literature review showed that recent studies prefer to apply ensemble models. Hong et al. (2018) applied ensemble models, namely AdaBoost, Bagging and Rotation Forest with J48 decision tree (JDT) as the base classifier, for LS in Jiangxi province of China. Their comparative analysis showed that ensemble models present higher accuracy in comparison to classical single model, i.e. JDT-RF had greater prediction capability that JDT or RF. In the same way, recently, Sachdeva et al. (2020) compared ensemble models comprising LR, Gradient Boosted Decision Trees and Voting Feature Interval for LSM in Brahmaputra valley region. They also found similar result that ensemble LR-GBDT-VFI offer higher accuracy that conventional techniques.
In the present analysis, the RFC model had the highest accuracy (AUC ¼ 97.0%). However, FDEMATEL-ANP (AUC ¼ 92.4%) and LR (AUC ¼ 91.2%) models outperformed the XGBoost (AUC ¼ 87.1%) and NBC (AUC ¼ 90.9%) machine learning methods whereas FR model had the lowest predictive accuracy (AUC ¼ 86.1%). In the study conducted by Cao et al. (2020), XGBoost performed with the highest accuracy whilst Achour and Pourghasemi (2020) found that the RFC reached the highest accuracy. Both studies analyzed and compared different machine learning methods to identify the highest accurate and performant model.
The novelty of any study is very important to attract the attention of readers which also influence the study interest and future research scope. In LS mapping, most of the recent studies preferred to use machine learning models, either single model or ensemble models or comparison amongst different models. Instead of using only machine learning, our study focussed on fuzzy MCDM, bivariate, and multivariate analysis along with three machine learning models in order to determine the accuracy of other models in comparison to machine learning models. Indeed, there are hundreds of recent studies that applied single or ensemble machine learning models. Based on models applied, it is very difficult to find any novelty amongst various studies, but each study is novel in term of geographical location and study findings because different outcomes have been found by applying same models in different geographical location (Achour and Pourghasemi 2020;Cao et al. 2020;Dang et al. 2020;Tien Bui et al. 2020;Ali et al. 2021). The main goal of this study was to present a comparative analysis amongst fuzzy MCDM, bivariate statistics, multivariate statistics and machine learning algorithms. This comparative analysis found a different result in comparison to other studies that the fuzzy DEMATEL-ANP as a MCDM model had higher accuracy than some machine learning models like XGBoost and NBC. The present analysis proved that fuzzy DEMATEL-ANP model had better accuracy than bivariate and multivariate models and even deep learning XGBoost model. Therefore, it is highly recommended to implement selected heuristic methods also along with the machine learning models, instead of directly applying only machine learning models or its ensemble form.
As a result, in the time of machine learning (ML) application in LS and spatial prediction ( Shirzadi et al. 2018;Vakhshoori et al. 2019;Achour and Pourghasemi 2020;Cao et al. 2020;Dang et al. 2020;Tien Bui et al. 2020), the purpose of the current study was to compare the performance and prediction accuracy of other techniques also along with machine learning algorithms. The results showed that although all of these models had high rate of prediction accuracy (AUC¼ >80%), but still some heuristic models, particularly fuzzy DEMATEL-ANP can offer better accuracy than selected machine learning techniques, particularly, NBC and XGBoost.
Lastly, it is also important to discuss the advantages and disadvantages of applied models because each model has its own pros and cons. In general, the model performance depends on different research areas and factors related with it ). In the present study, RFC performed better, whereas fuzzy DEMATEL-ANP had higher accuracy than machine learning models, namely XGBoost and NBC. The fuzzy DEMATEL-ANP is an extension form of classical set theory and analytic hierarchy process. It has many advantages like it takes insufficient information and evaluation into account, it allows imprecise information and solve problems with higher complexity, whereas it has also some disadvantages because sometimes it becomes difficult to develop fuzzy as it requires several simulations before being able to be used in the real world (Velasquez and Hester 2013). The FR model helps in deducing spatial relationship between locations of landslide occurrence and conditioning factors, which improves the accuracy of classification (Lee and Pradhan 2007), but it is a traditional method to evaluate LS. LR is a simple machine learning algorithm which provides greater training efficiency and lesser over-fitting problem in many cases, but it requires independent variables that are linearly related to the log odds. The RF has high tolerance in terms of algorithms, outliers and noises which can process multi-dimensional data without over-fitting, but the main advantage of RF is that a greater number of trees can make the algorithm slow and unproductive for real-time predictions ). The main advantage of using NBC is that it is suitable for solving multi-class prediction problems, which works very fast and saves time, whereas NBC assumes that all variables are independent which faces zero-frequency problem, so the prediction through NBC can be wrong in some cases. The XGBoost method supports linear as well as CART as its classifier and it can better prevent over-fitting that also reduce the complexity in model, but XGBoost lacks smoothness . Therefore, looking towards the advantages and performance of each model, the selected models were applied to make a comparative analysis amongst the MCDM, bivariate, multivariate and machine learning algorithms.

Conclusion
Landslides are the most hazardous disasters to have made their existence in the Kysuca river basin. There are numerous studies, which have been conducted worldwide to predict accurate location and impact of landslides. However, there remains an extensive scope to be thought that could offer a practical view related to LS with adequate information for future study. Thus, the main aim in this study was to develop and compare six heterogeneous models, particularly FDEMATEL-ANP, FR, LR, RFC, NBC and XGBoost for LSM. The main goal was to compare traditional models with machine learning to explore classification accuracy and success rate because most of the previous studies preferred to use statistical models whilst recent studies focussed on the use of only machine learning algorithms for LS. Such comparison of heterogeneous models is still unexplored in LS assessment. Initially, a total of 1000 landslide pixels were collected and randomly divided into 70/30% proportion for training and validation. In addition, 16 landslide conditioning factors were identified and selected for modelling. However, after multi-collinearity test, four insignificant variables were removed. The computed models were validated using the ROC-AUC curves, specificity, sensitivity, and overall accuracy. The study revealed that the RFC is the most promising machine learning model with 97% accuracy (ROC-AUC 0.97), 99.1% sensitivity rate and 97.5% specificity rate. Further, the present study also revealed that fuzzy DEMATEL-ANP outweighed bivariate frequency ratio, multivariate logistic regression and some machine learning models such as NBT and XGBoost. It is proved in the present analysis that always all machine learning models not performed well in susceptibility mapping and spatial prediction. So, it is recommended to test some heuristic models before applying machine learning models. The predicted models are promising for land-use planning and developing hazard mitigation strategies in the future. Since the RFC was the overall best method, it is suggested to be implemented for a regional level landslide hazard mitigation plan on the basis of the developed RFC LS map. This study makes available a competent method to LS models. Further heuristic techniques and deep learning methods in LS analysis will be explored in future research.