Landslide Detection and Susceptibility Modeling on Cameron Highlands (Malaysia): A Comparison between Random Forest, Logistic Regression and Logistic Model Tree Algorithms

: We used remote sensing techniques and machine learning to detect and map landslides, and landslide susceptibility in the Cameron Highlands, Malaysia. We located 152 landslides using a combination of interferometry synthetic aperture radar (InSAR), Google Earth (GE), and ﬁeld surveys. Of the total slide locations, 80% (122 landslides) were utilized for training the selected algorithms, and the remaining 20% (30 landslides) were applied for validation purposes. We employed 17 conditioning factors, including slope angle, aspect, elevation, curvature, proﬁle curvature, stream power index (SPI), topographic wetness index (TWI), lithology, soil type, land cover, normalized di ﬀ erence vegetation index (NDVI), distance to river, distance to fault, distance to road, river density, fault density, and road density, which were produced from satellite imageries, geological map, soil maps, and a digital elevation model (DEM). We used predictive value (NPV), sensitivity, positive predictive value (PPV), speciﬁcity, root-mean-squared error (RMSE), accuracy, and area under the receiver operating characteristic (ROC) curve (AUC). Our results indicated that the AUC was 92%, 90%, and 88% for the LMT, LR, and RF algorithms, respectively. To assess model performance, we also applied non-parametric statistical tests of Friedman and Wilcoxon, where the results revealed that there were no practical di ﬀ erences among the used models in the study area. While landslide mapping in tropical environment such as Cameron Highlands remains di ﬃ cult, the remote sensing (RS) along with machine learning techniques, such as the LMT model, show promise for landslide susceptibility mapping in the study area.

Due to difficult environmental conditions in the tropics, improved landslide detection and field surveys are critical for producing better LSMs. Herein, we firstly detected landslides over the Cameron Highlands in Malaysia using InSAR and then we prepared landslide susceptibility maps by three algorithms: LMT, RF, and LR. The LR is a functional algorithm based on statistical theory which has been used increasingly for spatial prediction of landslides in regions throughout the world. The LMT and RF models have recently been applied for this purpose and have different results in different case studies due to differences in probability distribution function (PDF). The main difference between this study and many other landslide studies is our application of both landslide detection using RS data, and landslide susceptibility mapping using strong machine learning algorithms, to prepare accurate maps of susceptibility to better manage landslide-prone areas in the study area. The RS analysis, modeling process, and final landslide maps were conducted in ENVI 8.00, WEKA 3.6.9 and ArcGIS 10.2 software, respectively.

Study Area
We selected an area of the Cameron Highlands, Malaysia, known for landslide activity for mapping and assessing landslide susceptibility. The study area covers an area of 81,250 km 2 and is located in the southwestern part of the Cameron Highlands. Most landslides here are attributed to land clearing for settlements and plantations ( Figure 1) [6]. The Cameron Highlands is a unique district in Pahang State. The Cameron Highlands is located in the eastern flank of the Malaysian geological main range, which is mainly composed of granite [79,80]. The granite in the study area is classified as megacrystic biotite granite [81]. In the scope of the study geological formations fall into two main groups: (1) acid intrusive, which covers most of the study area, and (2) the Ordovician-Silurian formation comprised of schist, phyllite, slate, and limestone. The highest and the lowest elevations of the study area are 1944 m and 953 m, respectively [82]. In tropical areas like the study area, average annual rainfall is quite high.  In the scope of the study geological formations fall into two main groups: (1) acid intrusive, which covers most of the study area, and (2) the Ordovician-Silurian formation comprised of schist, phyllite, slate, and limestone. The highest and the lowest elevations of the study area are 1944 m and 953 m, respectively [82]. In tropical areas like the study area, average annual rainfall is quite high. Based on the data of the Tropical Rainfall Measuring Mission (TRMM), the average annual rainfall recorded in the Cameron Highlands ranges between 3800 mm to 4200 mm. Malaysia has two main rainy seasons: February to May and September to December [79]. Between March and May and from November to December, rainfall in the study area reaches its peak [83]. Heavy rainfall flowing into the rivers during this time period resulted in landslides [79]. Intensity and amount of precipitation is another triggering factor that affects the slopes during such times.

Database and Data Collection
In preparing landslide conditioning factors, the construction of the spatial database and data collection, such as soil and geological maps, rainfall data, digital elevation model (DEM), and satellite imagery, are very important [6,84,85]. In the current study we used an AirSAR DEM (10 m spatial resolution) from Universiti Teknologi Malaysia (UTM), from which we produced seven different factors of slope, aspect, elevation, curvature, profile curvature, stream power index (SPI), and topographic wetness index (TWI) in ArcGIS. We used a geology map from Malaysias Mineral and Geoscience Department, from which the lithology and fault lines were digitized for the study area. From Malaysian Department of Agriculture, we acquired a soil map and digitized the soil layer. We downloaded the main road network from Open Street Map (OSM), while we generated rainfall data from TRMM data. We obtained Sentinel-2 satellite images from Sentinels Scientific Data Hub (https://scihub. copernicus.eu/) [86] and used it for generating the normalized difference vegetation index (NDVI) layer. Finally, we amassed Sentinel-1 and Landsat-8 satellite images from the Sentinels Scientific Data Hub (https://scihub.copernicus.eu/) website [86], and combined these to extract the land cover data. Table 1 highlights the factors we used in this study and their source and scale.

Research Methodology
Conditioning factors (e.g., slope) cannot directly cause landslide occurrences, while those which can create landslides (e.g., earthquakes) are triggering factors. The first step of the current research was preparation of the conditioning and triggering factors using GIS, System for Automated Geoscientific Analyses (SAGA), and ENvironment for Visualizing Image (ENVI) software. The conditioning and triggering factors were classified and then reclassified into several subclasses ranging from 2 to 9. The convergence of the InSAR, GE, and field survey techniques were used for landslide inventory by applying Sentinel Application Platform (SNAP) software, ENVI software, and GIS. A database of 152 landslides was mapped, of which 80% (122 landslide locations) were randomly selected for the training dataset, and the remaining 20% (30 landslide locations) were selected as the validation dataset. Three algorithms of LMT, LR, and RF were employed for application of landslide susceptibility mapping in the study area. The results were validated and compared by parametric and non-parametric models used in the study. Figure 2 represents the exact methodology used.

Research Methodology
Conditioning factors (e.g., slope) cannot directly cause landslide occurrences, while those which can create landslides (e.g., earthquakes) are triggering factors. The first step of the current research was preparation of the conditioning and triggering factors using GIS, System for Automated Geoscientific Analyses (SAGA), and ENvironment for Visualizing Image (ENVI) software. The conditioning and triggering factors were classified and then reclassified into several subclasses ranging from 2 to 9. The convergence of the InSAR, GE, and field survey techniques were used for landslide inventory by applying Sentinel Application Platform (SNAP) software, ENVI software, and GIS. A database of 152 landslides was mapped, of which 80% (122 landslide locations) were randomly selected for the training dataset, and the remaining 20% (30 landslide locations) were selected as the validation dataset. Three algorithms of LMT, LR, and RF were employed for application of landslide susceptibility mapping in the study area. The results were validated and compared by parametric and non-parametric models used in the study. Figure 2 represents the exact methodology used.

Landslide Inventory
Cost-effective information about landslide occurrence and processes can be extracted from satellite images [85]. For mapping areas susceptible to landslides, an accurate and reliable landslide detection system is very important [87]. It is generally accepted that the most exact way to identify and characterize landslides is through field investigation, which is a costly and time-consuming process [3]. Using satellite imagery, however, changes between two different images taken in separate positions and days can be recorded [86]. Two or more satellite images of the same scene and also the same position are needed to generate a land deformation map using InSAR technique [88,89]. From the http://scihub.copernicus.eu website [86], two imageries (Sentinel-1) of the different positions and dates were acquired for creating a map of interferogram. The technical characteristics

Landslide Inventory
Cost-effective information about landslide occurrence and processes can be extracted from satellite images [85]. For mapping areas susceptible to landslides, an accurate and reliable landslide detection system is very important [87]. It is generally accepted that the most exact way to identify and characterize landslides is through field investigation, which is a costly and time-consuming process [3]. Using satellite imagery, however, changes between two different images taken in separate positions and days can be recorded [86]. Two or more satellite images of the same scene and also the same position are needed to generate a land deformation map using InSAR technique [88,89]. From the http://scihub.copernicus.eu website [86], two imageries (Sentinel-1) of the different positions and dates were acquired for creating a map of interferogram. The technical characteristics of images used in  Table 2, and the Sentinel-1 images and location of the study site are shown in Figure 3. of images used in this study are listed in Table 2, and the Sentinel-1 images and location of the study site are shown in Figure 3.  InSAR utilizes the differences of phase values extracted from two different images [5,90]. Because of shorter wavelength of the C-band (compared with L-bands) the C-band SAR data cannot penetrate dense vegetation. Given that our study area has both densely vegetated and cleared areas, C-band SAR data is suitable for InSAR analysis. Before using InSAR, raw imageries need to be preprocessed. In this research, SNAP software was utilized to pre-process and process the imageries. The details of the methodology are illustrated in Figure 4. InSAR utilizes the differences of phase values extracted from two different images [5,90]. Because of shorter wavelength of the C-band (compared with L-bands) the C-band SAR data cannot penetrate dense vegetation. Given that our study area has both densely vegetated and cleared areas, C-band SAR data is suitable for InSAR analysis. Before using InSAR, raw imageries need to be pre-processed. In this research, SNAP software was utilized to pre-process and process the imageries. The details of the methodology are illustrated in Figure 4.
There are three sub-swaths in the Single Look Complex (SLC) product type of Sentinel-1 imageries; based on the study area, the suitable sub-swath should be split [91]. Once the appropriate sub-swath is selected, the orbit file must be updated. One of the primary processing tasks of InSAR is accurate co-registration of master and slave imageries. The slave imagery is then resampled to meet the geometry of the master one [92,93]. Once the imageries are well-stacked, the co-registered imageries need to be enhanced spectrally.
Cross multiplication of imagery pixels is used to form the interferogram map [92,93]. An interferogram output is created by measuring the phase difference between two SAR images [91,94] (Figure 5a,d). There are some bursts in Terrain Observation by Progressive Scans (TOPS) data, meaning that the imagery is divided into a few locations and any information within these locations must be corrected (Figure 5c). The curvature of the Earth results in some topographic distortion; therefore, the interferogram output must be flattened by removing the topographic phase [91,92]. A proper filter must be performed to remove noise from the interferogram [91]. The Goldstein filter produces the best noise reduction [5].  Step-by-step for assessment of landslides in study.
There are three sub-swaths in the Single Look Complex (SLC) product type of Sentinel-1 imageries; based on the study area, the suitable sub-swath should be split [91]. Once the appropriate sub-swath is selected, the orbit file must be updated. One of the primary processing tasks of InSAR is accurate co-registration of master and slave imageries. The slave imagery is then resampled to meet the geometry of the master one [92,93]. Once the imageries are well-stacked, the co-registered imageries need to be enhanced spectrally.
Cross multiplication of imagery pixels is used to form the interferogram map [92,93]. An interferogram output is created by measuring the phase difference between two SAR images [91,94] ( Figure 5A,D). There are some bursts in Terrain Observation by Progressive Scans (TOPS) data, meaning that the imagery is divided into a few locations and any information within these locations must be corrected ( Figure 5C). The curvature of the Earth results in some topographic distortion; therefore, the interferogram output must be flattened by removing the topographic phase [91,92]. A proper filter must be performed to remove noise from the interferogram [91]. The Goldstein filter produces the best noise reduction [5]. Step-by-step for assessment of landslides in study.

Figure 4.
Step-by-step for assessment of landslides in study.
There are three sub-swaths in the Single Look Complex (SLC) product type of Sentinel-1 imageries; based on the study area, the suitable sub-swath should be split [91]. Once the appropriate sub-swath is selected, the orbit file must be updated. One of the primary processing tasks of InSAR is accurate co-registration of master and slave imageries. The slave imagery is then resampled to meet the geometry of the master one [92,93]. Once the imageries are well-stacked, the co-registered imageries need to be enhanced spectrally.
Cross multiplication of imagery pixels is used to form the interferogram map [92,93]. An interferogram output is created by measuring the phase difference between two SAR images [91,94] ( Figure 5A,D). There are some bursts in Terrain Observation by Progressive Scans (TOPS) data, meaning that the imagery is divided into a few locations and any information within these locations must be corrected ( Figure 5C). The curvature of the Earth results in some topographic distortion; therefore, the interferogram output must be flattened by removing the topographic phase [91,92]. A proper filter must be performed to remove noise from the interferogram [91]. The Goldstein filter produces the best noise reduction [5]. By increasing the integer number of cycles, the phase difference can be reduced [92]. Herein, 2π is performed to separate the range of the phase unwrapping [5,90]. Whenever time in unwrapping process reaches 2π, the cycle should be repeated [90,92]. There is always ambiguity in 2π; to address it, we applied the phase unwrapping process using SNAPHU Software (Figure 5b). One of the standard processes in the InSAR technique is terrain correction, by which the generated interferogram is corrected geometrically [92,93].

Landslides Detection by Overlying the Phase, Unwrapped and Coherence Bands
Google Earth (GE) was first introduced in 2001, with outstanding advantages using image data obtained from satellites and aircraft (i.e., remote sensing technology) [95,96]. By zooming into GE, the finer details appear [87]. Before the advent of Landsat-8 satellite data, the images of Landsat-7 were employed to cover most of the GE images. The spatial resolution of the GE images is 15 m [95,96]. Satellite images with global spatial coordinates are captured and shown on GE [95,97]. In the present study, we applied the phase, unwrapped, and coherence bands, with the assist of GE, to identify and record historical landslide locations. We used GE data to help identify landslides and then enlarged satellite images using the phase, unwrapped, and coherence bands, the landslide locations were identified and digitized. The soil and geology in the study area tend to host translational debris slides and creeping failures.

Validation of the Detected Landslide Locations
Generally, an assessment of accuracy denotes the quality of the obtained information [98]. Global Positioning System (GPS) technology is used in many deformation studies to identify locations of features [99,100]. In the current study, we used a handheld GPS to validate the extracted landslide locations. The type of antenna will affect the spatial-resolution of a GPS [99,101]. Because of the hilly terrain of the study area, we field checked 20% of the detected landslide locations, near the main road, for accuracy.

Topographical and Geomorphological Factors
DEM is a digital representation of the earth surfaces [3,85,102,103]. A 10-m DEM was directly used to create a slope, aspect, elevation, curvature, profile curvature, SPI and TWI using ArcGIS, and SAGA software. Additionally, by using the hydrology toolbox in ArcGIS the stream networks of the study area were delineated from DEM and were utilized to produce the proximity to river and the river density layers.

Lithology
Lithology determines the degree of strength and permeability of soils and rocks; therefore, geology greatly influences slope stability [25,[104][105][106]. Due to high infiltration rate and reduced cohesion, unconsolidated formations are susceptible, with less shear strength than consolidated ones [11,28,107,108]. Further, surfaces near active fault lines tend to be fractured, allowing percolation of water and increased rates of weathering. Faults tend to reduce rock strength and can predispose slopes to landslides [3,109,110]. In this study the layers of lithology and fault lines were digitized from a 1:100000 scale geologic map. The acid intrusive formation (61 km 2 out of 81.250 km 2 ) and the Ordovician-Silurian formation are the only two formations in the study area.

NDVI
Evapotranspiration, rainfall, and infiltration have a remarkable impact on vegetation coverage, which itself has substantial effects on the soil characteristics [105,106,111,112]. The incidence of landslides is reduced by the increased soil strength and root reinforcement afforded by vegetation [11,19]. The NDVI map for the present study was generated from Sentinel-2 satellite data acquired on 11 October 2017, which was calculated using ArcGIS by the common formula of Float (NIR-Red)/ (NIR + Red).

Land Cover
Many works of literature have focused on the impact of land cover on slope stability [6,22,27,[113][114][115][116][117][118]. Vegetation coverage increases the stability and durability of slopes by providing shear strength, absorbing soil moisture, and reducing erosion [18,[119][120][121]. Surfaces with dense vegetation are less prone to landslides compared to sparse surfaces [6,105,106]. In the current study, the land cover map was extracted from the combination of the Landsat-8 and the Sentinel-1 satellite imageries. Five land covers of florification, forest, water body, cleared forest, and township were mapped and used.

Road Networks
Road networks are anthropologically induced instabilities, and can be an important source of landslides in the mountainous regions [10,104,122,123]. OSM was used to download the main road networks of the study area. The entire road network in the study area, about 32 km in length, was applied to create the proximity to road and road density layers.

Soil Type
The influences of soil on landslides have been recognized in previous studies, and researchers have widely used soil layer as one of the main parameters to locate areas susceptible to landslides [29,119,124,125]. The difference between shallow and deep landslides depends greatly on the thickness and material of soil [85,120,126,127]. In this study, the soil map was obtained from the Department of Agriculture, Malaysia and used for digitizing the soil factor. There are only two soil types in the study area: the serong series and the alluvium-colluvium type.

Rainfall
Malaysia is a tropical country, with heavy rainfall almost year-round. TRMM data was used to create the rainfall layer of the study area. Monthly data for all months of 2017 was obtained and used to create average annual data through ArcGIS. The rainfall recorded in 2017 ranged from 3800 mm to 4200 mm.
The classification of environmental factors is one of the main steps in any landslide susceptibility assessment [6,31,32,105,128]. Taking into account the mechanism of the landslide occurrence, previous studies (local or global), and the characteristics of the study area, 17 landslide conditioning and triggering factors were identified and classified ( Table 3). The conditioning factors used in this study is shown in Figure S1 as supplementary materials. These classifications are associated with the evaluation of each factor's sensitivity to landslide occurrence [6,85,105,106]. Table 3. Classes of environmental factors affecting landslide susceptibility.

Logistic Regression (LR) Model
The LR model has been widely utilized to analyze binary variables [125,129,130]. It describes the independent parameters and the dependent parameter, where the relationship between the dependent variable and the independent variables is nonlinear. In this study, the relationship between landslide occurrence and the landslide-affecting factors was described by the LR model, which Equations (1) and (2) express well: where, landslide occurrence's probability is described by p (0 ≤ p ≤ 1), at the same time, Z defines linear logistic parameter (−∞ ≤ p ≤ +∞). Equation (2) further analyzes the probability: where, n represents the quantity of landslide factors and b 0 is the constant coefficient.

Logistic Model Tree (LMT)
Because it combines the decision tree and linear logistic regression methods, the Logistic Model Tree (LMT) is a comprehensive algorithm [131,132]. In comparison with a traditional decision tree, the logistic regression function is used by LMT to value the probability of the classes. It utilizes the LogitBoost for creating the logistic regression functions [31,131]. Equation where Hc (x), the least-squares fit, is transformed so that c c=1 Hc (x) = 0, and c is the number of landslides in each class.

Random Forest (RF) Model
The RF model is well-known as a robust AI technique for classification of natural hazards, especially landslides [31,121]. It was introduced first by Tin Kam Ho [133] and was developed by Breiman in 2001 [134]. In the RF model, decision trees are developed based on the bootstrap samples of the dataset used [135]. For each sample, only one decision tree is built. In each tree, the nodes are defined and transferred to the leaves [136].
In this study, the RF consists of two trees (landslide and non-landslide), each constructed using random features. Typically, in a RF algorithm, the generalization error is expressed as follows [137]: where x and y are landslide conditioning factors that indicate the probability over x and y space, mg is the margin function, and I( * ) is the indicator function [134,137]. Advantages of the random forest model include low computational cost, good performance in large data, utilization of a large amount of input features without feature elimination, determination of the most important parameters in classification, identification of relationships between factors and classification, and avoidance of over-fitting [138].

Evaluation Methods
One of the essential purposes of comparing models is to evaluate their performance, which determines whether a model is applicable or not [6,105].

Statistical Measurements
In this research, the ROC curve and statistical measurements (i.e., specificity, sensitivity, accuracy, PPV, NPV, and RMSE), and the Friedman and Wilcoxon tests were applied to assess the models' performances [139,140]. Table 4 shows the confusion matrix of four types of possible landslides: true positive (TP), false positive (FP), false negative (FN), and true negative (TN). The TP denotes the landslides points, while TN denotes the non-landslides points. The FN and FP are the numbers of pixels that incorrectly classified as non-landslides and landslides, respectively [141]. The sensitivity indicates the number of landslides that are accurately classified on the total number of predictions (Equation (5)). Specificity refers to the ratio of the number of incorrectly classified cells to the total predicted non-landslide cells (Equation (6)). Accuracy is the rate of landslide and non-landslide pixels which were correctly classified (Equation (7)) [6,142]. The error between the predicted and observed data is indicated by RMSE (Equation (8)) [6,29,142].
Speci f icity = TN TN + FP where the total training dataset defined by N, X predicted and X actual are the predicted sample in the training data and the actual value from the landslide algorithms, respectively.

Receiver Operating Characteristics (ROC) Curve
The ROC curve is an applicable way to evaluate the validity and reliability of different models used in a variety of studies, including landslide assessment [11,102,111,143]. The ROC curve and the area under the ROC curve (AUC) are widely employed as effective methods for validation of landslide studies [6]. The ROC curve plots the true positive rate, correctly classified landslides, or sensitivity, against true negative rates, correctly classified non-landslides or specificity [6,144]. The results of the ROC curve are concluded through the AUC, in which a classification implies as AUC = 0.5, and a complete success represents as AUC = 1 [6,145,146]. The statistical measurements of AUC are as Equation (9) [147]; where N and P denote the number of non-landslides and landslides points, respectively.

Friedman and Wilcoxon
Demšar [148] noted that, because non-parametric tests such as the Friedman and the Wilcoxon do not assume homogeneity of variance or normal distribution, their results are more robust than parametric tests [6,149,150]. In the present study, the Friedman and Wilcox tests were employed to validate the performance of the algorithms used. The Friedman test is based on the two assumptions of H 0 and H 1 (p-values) [6,131,149,150]. In the first assumption (H 0 ), there is no significant differences between the models applied, while in the second assumption there is clear statistical difference between the models. The Friedman test is based on group comparison rather than pairwise [6]. Since the Wilcoxon signed-rank model considers pairwise comparison for every single model, we can clearly see the performance of every model compared to the another. In this model, the p-value (significance) and z-value are utilized for assessing the rate of differences among the algorithms [6,132]. A p-value less than 0.05, and a z-value greater than −1.96 and +1.96, indicate the performance of the susceptibility algorithms is different [6,131,132].

Landslide Inventory
A fast and accurate way to detect historical landslides [151] is to use the remote sensing (RS) technique [152]. In this study, the historical landslides were identified using InSAR remote sensed imagery, and overlaying the phase, unwrapped, and coherence bands on the Google Earth (GE) images. Detected landslide locations on the different phase, unwrapped, and coherence bands are shown in Figure 6. Note that we used Sentinel-1 satellite imageries and InSAR based on SNAP software. In the phase and unwrapped bands, the red color shows the highest values, whereas the blue color represents the lowest values (Figure 6a,b, respectively). White surfaces (high values) are distinguished from stable areas (dark values), and deformed surfaces are shown better in the combined band (Figure 6c).
In the inventory stage, we extracted the phase, unwrapped and coherence groups from the InSAR map. Subsequently, we converted them into the KML/KMZ format. We identified and recorded landslides by using surface deformation maps (phase, unwrapped, and coherence bands) and zooming in and scrolling around the deformed areas. By zooming in on the GE images, the finer details appear; therefore, the deformed surfaces become more visible (especially in vegetated areas, where landslides left exposed soil). Since much of the study area is sparsely vegetated, this method is suitable for landslide detection (Figure 7).
Forests 2020, 11, x FOR PEER REVIEW 13 of 28 blue color represents the lowest values (Figure 6a and 6b, respectively). White surfaces (high values) are distinguished from stable areas (dark values), and deformed surfaces are shown better in the combined band (Figure 6c). In the inventory stage, we extracted the phase, unwrapped and coherence groups from the InSAR map. Subsequently, we converted them into the KML/KMZ format. We identified and recorded landslides by using surface deformation maps (phase, unwrapped, and coherence bands) and zooming in and scrolling around the deformed areas. By zooming in on the GE images, the finer details appear; therefore, the deformed surfaces become more visible (especially in vegetated areas, where landslides left exposed soil). Since much of the study area is sparsely vegetated, this method is suitable for landslide detection (Figure 7).

Validation
We used a handheld GPS to validate 20% (30 locations) of the detected landslide locations in the field, to verify whether the landslides were correctly recorded. Figure 8 shows the geographical position of the validated landslide locations, 64 landslide locations out of the 152 features have an area larger than 500 m 2 , and the remaining 88 locations are smaller than this value. Overall, the

Validation
We used a handheld GPS to validate 20% (30 locations) of the detected landslide locations in the field, to verify whether the landslides were correctly recorded. Figure 8

Validation
We used a handheld GPS to validate 20% (30 locations) of the detected landslide locations in the field, to verify whether the landslides were correctly recorded. Figure 8 shows the geographical position of the validated landslide locations, 64 landslide locations out of the 152 features have an area larger than 500 m 2 , and the remaining 88 locations are smaller than this value. Overall, the smallest and the largest landslides have areas of 40 m 2 and 5735 m 2 , respectively.

Generating Landslide Susceptibility Mapping (LSM)
In this study we applied three algorithms of LR, LMT, and RF to the landslide susceptibility. We used 17 conditioning variables for training the models. We measured the landslide susceptibility indices by multiplying the weights (calculated through Waikato Environment for Knowledge Analysis [WEKA] software) by the secondary reclassified conditioning factors. Based on our findings, all the models yielded promising algorithms for LSM in the study area. Figure 9 shows the LSM generated by the LR model, which we ultimately reclassified into four susceptibility classes of low, moderate, high, and very high using the natural break method in ArcGIS. According to the LR susceptibility map, the very high susceptible zones comprised approximately 22% (18.2 km 2 ) of the total area (81.2 km 2 ), while roughly 23% (18.5 km 2 ) was categorized as high susceptibility and 23% of the total area (18.9km 2 ) was classified as moderate susceptible zone. Finally, 32% (25.7 km 2 ) belonged to the low susceptibility class. Figure 9 shows the LSM generated by the LR model, which we ultimately reclassified into four susceptibility classes of low, moderate, high, and very high using the natural break method in ArcGIS. According to the LR susceptibility map, the very high susceptible zones comprised approximately 22% (18.2 km 2 ) of the total area (81.2 km 2 ), while roughly 23% (18.5 km 2 ) was categorized as high susceptibility and 23% of the total area (18.9km 2 ) was classified as moderate susceptible zone. Finally, 32% (25.7 km 2 ) belonged to the low susceptibility class.

LSM using the Logistic Model Tree (LMT) Model
The LSM produced from the LMT model is shown in Figure 10. The final extracted susceptibility map from this model was reclassified into four different susceptibility zones. The results show that the moderate and high susceptibility zones yielded about 23% (roughly 18 km 2 ) of the total area each. The low susceptible zone comprised 30% (24.7 km 2 ), while for the very high susceptible areas this value was 24% (19.1 km 2 ).

LSM using the Logistic Model Tree (LMT) Model
The LSM produced from the LMT model is shown in Figure 10. The final extracted susceptibility map from this model was reclassified into four different susceptibility zones. The results show that the moderate and high susceptibility zones yielded about 23% (roughly 18 km 2 ) of the total area each. The low susceptible zone comprised 30% (24.7 km 2 ), while for the very high susceptible areas this value was 24% (19.1 km 2 ). The RF algorithm fell into four susceptibility zones of low, moderate, high, and very high ( Figure  11). However, based on the findings, 27% (22.0 km 2 ) and 26% (21.2 km 2 ) of the total area fell into the

LSM by the Random Forest (RF) Model
The RF algorithm fell into four susceptibility zones of low, moderate, high, and very high ( Figure 11). However, based on the findings, 27% (22.0 km 2 ) and 26% (21.2 km 2 ) of the total area fell into the high and very high susceptibility zones, respectively. These values were slightly lower for the low susceptibility zone value of 25% (20.0 km 2 ). Finally, the moderate areas yielded about 22% (18.1 km 2 ) from the total area of 81.2 km 2 .

LSM by the Random Forest (RF) Model
The RF algorithm fell into four susceptibility zones of low, moderate, high, and very high ( Figure  11). However, based on the findings, 27% (22.0 km 2 ) and 26% (21.2 km 2 ) of the total area fell into the high and very high susceptibility zones, respectively. These values were slightly lower for the low susceptibility zone value of 25% (20.0 km 2 ). Finally, the moderate areas yielded about 22% (18.1 km 2 ) from the total area of 81.2 km 2 . Figure 11. Landslide susceptibility map produced from the random forest (RF) model.

Statistical Measurements
Performance results of the models from the validation dataset are shown in Table 5. The validation phase shows that the highest values of PPV and NPV belong to the LMT model, followed by the LR and RF models. Furthermore, results indicated that the LMT model had the highest sensitivity (68.75%), meaning 68.75% of the landslides pixels were correctly classified in the landslide classes, followed by the LR and RF models at 65.63% and 62.50%, respectively. The highest specificity (71.43%) belonged to the RF model, indicating 71.43% of the non-landslide pixels were correctly classified, followed by the LR and LMT models at 67.86% and 64.29%, respectively. In comparison with the LR model, the RF and LMT models had the highest accuracy value of 70% and 66.67%, respectively. The RMSE value was 0.3, 0.4 and 0.4 for the LMT, LR and RF models, in order. Overall, all models successfully trained with the validation datasets; therefore, the models can extract maps of areas susceptible to landslides in the study area.

ROC Curve
The capability of the prediction accuracy of the models used was also evaluated using the AUC and validation dataset. The comparison of AUC for the models is shown in Figure 12. The prediction rate curve calculated the prediction capability of the models [6]. The higher AUC values of 0.92 and 0.90, were obtained for the LMT and the LR algorithms, respectively, while an AUC value of 0.88 was calculated for the RF model. Overall, it can be concluded that all models had a high prediction rate, with few differences between them for mapping areas prone to landslides in the study area.

ROC Curve
The capability of the prediction accuracy of the models used was also evaluated using the AUC and validation dataset. The comparison of AUC for the models is shown in Figure 12. The prediction rate curve calculated the prediction capability of the models [6]. The higher AUC values of 0.92 and 0.90, were obtained for the LMT and the LR algorithms, respectively, while an AUC value of 0.88 was calculated for the RF model. Overall, it can be concluded that all models had a high prediction rate, with few differences between them for mapping areas prone to landslides in the study area.

The Friedman and Wilcoxon
To evaluate the algorithms, two non-parametric statistical methods of the Friedman and Wilcoxon signed-rank were also applied. Table 6 highlights the result of the Friedman test. It revealed the mean rank values for the LMT, LR, and RF models of 2.07, 1.73, and 2.20, respectively. Results

The Friedman and Wilcoxon
To evaluate the algorithms, two non-parametric statistical methods of the Friedman and Wilcoxon signed-rank were also applied. Table 6 highlights the result of the Friedman test. It revealed the mean rank values for the LMT, LR, and RF models of 2.07, 1.73, and 2.20, respectively. Results indicated significant (Sig.) differences among the models because of significance lower than 0.05 (Sig. = 0.030). However, the tests did not indicate between which models there were statistical differences. The Wilcoxon signed-rank test was employed to assess pairwise differences among the models (at 5% significance level). Once there is no remarkable difference among the models at 5% significance (rejection of the null hypothesis), it will be proved that the results are not the same [6,117]. Tien Bui et al. [117] defined that when p (value) is lower than 0.05 and z (value) higher than −1.96 and +1.96, the capability of the algorithms is expected to be significantly different. According to Table 7, the results of the Wilcoxon test concluded that there were no high statistical differences among the used algorithms.

Discussion
Landslide modeling and landslide susceptibility mapping (LSM) face challenges due to uncertainties in the quality and resolution of input data. Landslide researchers strive to improve the detection and mapping of landslide locations and conditioning factors [33]. We conducted this study in a tropical environment of Cameron Highlands, Malaysia. It was carried out in two parts: (1) landslide locations were detected by RS data and (2) landslide susceptibility modeling was done by three machine learning algorithms, namely LMT, LR, and RF. We selected 17 conditioning factors for landslide susceptibility modeling. Detection and recognition directly on the ground is generally a difficult task in conditions of dense vegetation, cloudy and rainy weather, and inaccessible mountainous terrain [6]. Improved detection of landslides in such environments can be performed using GIS and RS techniques. Cheng et al. [153] and Metternicht et al. [154] believe that these techniques can be useful for landside inventory mapping, landslide susceptibility, and hazard mapping. We effectively detected and mapped landslide locations using RS data, and they were validated by ground control points. Our findings indicated reasonable concordance between landslides observed and estimated over the study area, but we acknowledge that we may have missed some small landslides within densely forested areas as also found by Brardinoni et al. [155].
After checking the accuracy of landslides detected, we used three machine learning algorithms for landslide modeling and susceptibility assessment. Our findings showed that among the three landslide models, the LMT decision tree model outperformed and outclassed the LR and RF models. What is remarkable is the difference between modeling results due to the difference in structure and probability distribution function (PDF) of algorithms. There was a difference when two or more algorithms were applied and compared from one study area to another because of uncertainties in the data input and changes in geo-environmental factors. For example, although LR and RF are strong algorithms that in some studies were more powerful and applicable for prediction of landslides, in this study the LMT was more flexible and powerful. The LMT better decreased the noise and over-fitting problems during the tree pruning process by the CART algorithm, and hence, the performance and goodness-of-fit were higher than the other models. Tien Bui et al. [131] reported the LMT could be useful as a promising technique in shallow landslide modeling.
Additionally, Chen et al. [132] concluded that although the RF had the highest prediction accuracy, the LMT outperformed the CART algorithms for shallow landslide mapping. Moreover, the LR outperformed the RF decision tree model in the area. The LR model is a benchmark and statistical-based theory that had good results for spatial prediction of landslides in some studies [156][157][158][159].
Although RF showed the lowest prediction accuracy in comparison to the LMT and LR, the RF with AUC equal to 0.88 yielded an acceptable result. The RF is a powerful algorithm that in some studies showed more power than the other models. For example, Chen at al. [70] compared the efficiency of best-first decision tree, random forest, and naïve Bayes tree for landslide susceptibility modeling. They found that the RF had the highest AUC values (0.869), smallest standard error (0.025), narrowest 95% confidence interval (0.819-0.918), highest accuracy value (0.774), highest precision (0.662), and highest F-measure (0.662) for the training dataset. The RF algorithm, when more trees are grown and added during the modeling process by the training dataset, does not create over-fitting, but a small error is generated [160,161]. Zhang et al. [162] state that the RF produces high prediction accuracy even when using a large number of input variables without variable deletion, and returns very small sets of classifications. Overall, all three landslide models have given reasonable performance for landslide susceptibility assessment but, comparatively, the LMT model has given the best performance. Therefore, we propose the LMT as a promising technique for tropical environment in areas with similar attributes as the Cameron Highland for LSM.

Conclusions
Remote sensing and geographic information system technologies have been widely used in the application of different models to generate landslide susceptibility maps, which are very important for policy and decision-makers to monitor and protect areas prone to landslides. The current study released the results of comprehensive research on landslide susceptibility mapping in the Cameron Highlands, Malaysia. It is generally accepted that landslides are the most frequent and dangerous natural disasters in the Cameron Highlands. Due to the increase in urbanization and plantations, landslide occurrence is expected to rise in the study area. Authorities are concerned with the damages to property and human lives caused by natural disasters, including landslides. To better address the influences of natural disasters, policy-and decision makers have to manage the problems using a better understanding of landslide potential.
We used InSAR, GE imagery, and extensive field surveys to inventory landslides. In all, we identified and extracted 152 landslide and employed 80% (122 landslide locations) for training and the remaining 20% (30 landslide locations) for validation of the LR, LMT, and RF algorithms. We obtained and adapted the environmental factors from numerous databases, including a geological map, digital elevation model, soil map, satellite imagery, TRMM data, and OSM. From abovementioned databases we extracted different variables, including proximity to river, road and fault, road and river density, normalized difference vegetation index, elevation, slope, aspect, profile curvature, curvature, lithology, rainfall, land cover, soil, stream power index, and topographic wetness index. We validated the results using the ROC curve, statistical measurements, and the Friedman and Wilcoxon tests. Findings of the validation revealed there are few differences among the models for mapping areas susceptible to landslides in the study area. Using a GPS, 20% of the detected landslide locations were validated, and the results showed that the locations were correctly recorded. Results indicated the superiority of the LMT to LSM in the tropical environment, such as Cameron Highlands. We propose that in tropical environments where landslide mapping is difficult, RS data along with machine learning techniques, such as LMT model can be used as promising methods for LSM.