GIS-based evolution and comparisons of landslide susceptibility mapping of the East Sikkim Himalaya

ABSTRACT The main sought of this study is to assess the landslide susceptibility map (LSM) of the East Sikkim Himalaya based on the comparative model analysis using frequency ratio (FR), logistic regression (LR), random forest (RF) and integration of analytical hierarchy process (AHP) with FR (AHP-FR). The models were trained by 166 landslides (70% training) and 12 landslide causative factors whilst tested with the help of 71 landslides (30% testing). Their spatial correlation between the landslides and the causative factors was analysed by using a multicollinearity test. The generated LSM was classified into five classes, i.e. very low, low, moderate, high and very high. In East Sikkim, very high classes of the AHP-FR, LR and FR models cover the area of 11.97%, 11.99% and 7.13%, respectively. The accuracy of prepared LSM was evaluated by using the success rate curve (SRC), prediction rate curve (PRC) and seed calculation area index (SCAI). The area under the curve (AUC) of the success rate curve is 0.88 for the RF model, 0.85 for AHP-FR, 0.78 for LR and 0.79 for FR, whilst the prediction rate curve is 0.86% for the RF model, 0.81 for AHP-FR, 0.79 for LR and 0.78 for FR. The SCAI values of very high susceptibility classes are 0.14, 0.17, 0.18 and 0.19 for the RF, AHP-FR, LR and FR models, respectively. The RF and integrated AHP-FR methods indicate better results as compared to other statistical models.


Introduction
Landslide is often termed universally as a wide variety of processes and mass movements in a downward direction consisting of rock, earth and debris under the influence of gravity (Ayalew and Yamagishi 2005). It frequently occurs in hilly or mountainous regions and is extremely devastating to the environment and its inhabitants. It occurs when the shear stress of the material increases and the materials shear strength decreases on a slope. In 2005, World Bank prepared a report on global risk analysis, which reveals that approximately 3.7 million km 2 of land is susceptible to landslide, which could cause a risk of 300 million human lives (CRED 2019). India alone witnessed over 11,000 deaths due to landslides in the last 12 years and has also perceived the fastest rise in human-triggered fatal landslides globally Ray 2021). NASA reported about 10,804 rainfall-triggered landslides between 2007 and 2017 (Froude and Petley 2018;Kirschbaum and Stanley 2018;Petley 2012), whilst earthquake-induced landslides caused 55,997 fatalities (NASA 2020). According to Brabb (1993), 90% of landslide losses can be avoided if the problem has been recognized before the landslide event. Hence, there is a need for landslide susceptibility assessment at various spatial scales to forecast the phenomenon. The landslide susceptibility map is the degree of instability of the terrain that will show the likelihood of occurrences of future landslides, taking into consideration a set of influencing factors such as the slope/aspect of the terrain, the geology and geomorphology of the region (Basu and Pal 2020;Corominas et al. 2014); i.e. susceptibility = ƒ (Landslide, Landslide related factors). It predicts 'where' landslides are likely to occur (Guzzetti et al. 2006). Preparation of landslide susceptibility maps is the first step in managing the risk assessment in landslide-prone areas. It is the function of landslides and their causative factors.
There are various traditional methods for landslide susceptibility mapping based on the field visit and manual analysis. In recent years, the advancement of remote sensing techniques has helped to classify different triggering factors of landslides and identify the past landslides to prepare a landslide inventory map. It became very easy to use data in the GIS and produce landslide susceptibility mapping. However, GIS and remote sensing have become more popular in spatial analysis for the landslide (Zhao and Lu 2018;Anbalagan et al. 2015;Yalcin et al. 2011).
From the last few decades, many researchers have applied various methods for mapping of the LSM using the GIS technology, based on the qualitative and quantitative approaches (Corominas et al. 2014;Dai et al. 2002;Fell et al. 2005). The qualitative method includes a heuristic or index-based method, in which weights are assigned to the factors based on expert knowledge; hence, it is a subjective approach (Sarkar et al. 2008;Castellanos and Van Westen 2008;Sarkar and Kanungo 2004;Guzzetti et al. 1999;Pardeshi et al. 2013;Van Westen and Getahun 2003). The analytical hierarchy process and weighted linear combination are the two widely used heuristic models (Chen et al. 2016;Ganesh and Reddy 2014;Mishra and Sarkar 2020;Pourghasemi et al. 2013;Akgun et al. 2008;Ayalew et al. 2004). The quantitative methods based on the mathematical relationship between slope stability and causative factors are analysed to overcome the limitation of the heuristic method. The quantitative approaches also include geotechnical engineering (deterministic) analysis and statistical analysis. The deterministic method based on numerical models requires the detailed study of soil and different characteristics of slopes (Youssef et al. 2016). Because of the difficulty of acquiring data, this model is not accurate for landslide mapping at a regional scale (Yilmaz 2009). The statistical analysis is based on the numerical method that shows the relationship between causative factors and landslides. Bivariate and multivariate methods are the two different statistical methods (Pardeshi et al. 2013). Statistical methods include the frequency ratio (FR) model ; Lee and Pradhan 2007;Mondal and Mandal 2017), logistic regression (LR) method (Saro Lee 2005;Lombardo et al. 2015;Mishra and Sarkar 2020;Wang et al. 2019), weight of evidence (Arabameri et al. 2018;Chen et al., 2018a;Khosravi et al. 2016;Kumar and Anbalagan 2019;Ozdemir and Altural 2013;Pourghasemi et al. 2013), multivariate adaptive regression spline (Felicismio et al., 2013), multivariate regression (Akgun and Turk, 2011;Guzzetti et al. 2006), evidential belief function (Althuwaynee et al. 2015;Pradhan and Kim 2017), Certainty factor (Chen et al., 2018a;Devkota et al. 2013;Liu and Miao 2018;Wang et al. 2015), fuzzy logic (Kanungo et al. 2006;Champatiray et al. 2007) and information ratio model (Lin and Tung 2004;Saha et al. 2005). In recent times, various machine learning algorithms such as support vector machine (SVM) (Chang et al. 2019;Goetz et al. 2015;Wang et al. 2019;Zhang et al. 2018), artificial neural network (ANN) (Arora et al. 2004;HembRam et al. 2020;Wang et al. 2019), Naive Bayes Pham et al. 2017;Shirzadi et al. 2017), random forest (RF) (Ada and San 2018;Chang et al. 2019;Tang et al. 2020), decision tree (Pradhan 2013;Tien Bui et al. 2016;Vakhshoori et al. 2019) boosted regression trees (Youssef et al. 2015) have been applied in various fields of natural hazards. With the advancement of the machine learning approaches, deep learning has also been employed in different spatial mapping, viz., gully erosion (Chen et al. 2021b;Amiri et al. 2019), flood mapping (Lei et al. 2021;Shahabi et al. 2020;Pham et al. 2020b), piping erosion (Chen et al. 2021a) and groundwater potential (Chen et al. 2021b;Naghibi et al. 2017).
The literature review suggests that many methods have been developed for landslide susceptibility mapping; still, there is limited study regarding which type of model is appropriate for landslide mapping (Xu et al. 2012;Hong et al. 2016). Many studies are indicating the comparisons of heuristic (AHP) with different statistical methods. Pourghasemi et al. (2012) compared the fuzzy logic and AHP at Haraz watershed, Iran; it reveals that the fuzzy logic model performed better than AHP. Meena et al. (2019) showed a comparative study of FR, AHP and SMCE, in which FR and SMCE offer better results than AHP. Pourghasemi et al. (2013a) prepared LSM using binary logistic regression (BLR), AHP and statistical index (SI), their study indicates that the BLR and SI performed better than AHP. The above literature proved that the AHP model does not always show promising results for landslide susceptibility mapping. Therefore, it is necessary to improve the accuracy of AHP for susceptibility mapping. Recently, different authors attempted to integrate different models to improve the accuracy of landslides (Bui et al. 2016;Dou et al. 2020;Mondal et al., 2013). The present study's objectives are to improve the AHP model by integrating with FR and show the comparative analysis of the LR, FR, RF and AHP-FR for landslide susceptibility mapping in the East Sikkim area.

Study area
Sikkim Himalaya is a small mountainous state in the eastern Himalaya region spread in a 7096 km 2 area, bound by 27° 04ʹ46"N to 28°07ʹ48"N latitudes and 88°58"E to 88°55ʹ25" E longitudes with a maximum altitude of 8585 m. Due to tropical rain, rugged weather conditions, complex geology, thick soil cover and active tectonic, Sikkim Himalaya witnesses frequent landslides. Sikkim state comprises East Sikkim, North Sikkim, South Sikkim and West Sikkim. The present study is focused on the East Sikkim, which encompasses latitude 27°10ʹN -27°23ʹN and longitude 88°25ʹE -88°45ʹE ( Figure 1) and covers an area of about 942 km 2 . Nathula Hill is one of the main tourist attractions in East Sikkim, along with Tsongo Lake. Five major rocks of the Precambrian age cover the major portion of East Sikkim, i.e. Everest Pelitic formation, Rangoli schist formation, Darjeeling Gneiss Formation, Chungatang formation and Kanchanjunga gneiss. The area is traversed by two major thrusts, i.e. Chungatang thrust and Main Central Thrust. Chungthang formation includes quartz-biotite schist, Alsilicate rocks and graphite schists. The Kanchenjunga group comprises hard rocks, highly jointed and intruded by pegmatite and tourmaline granite at places.
The east Sikkim area experiences frequent natural and anthropogenic landslides, with a significant loss of properties and lives (Kaur et al. 2019). In 1968, hundreds of rainfall-induced landslides occurred in Sikkim and Darjeeling, resulting in the loss of 33,000 human lives, property and infrastructure. In 2011, the Sikkim earthquake triggered several landslides causing several casualties and properties (Singh et al., 2016). In the present study area, a total of 237 major landslide events were identified.

Materials and methods
Different thematic layers representing the slope, aspect, elevation, curvature, topographic witness index (TWI), distance from the road, distance from drainage, distance from lineament, rainfall, geology, land use land cover (LULC) and normalized difference vegetation index (NDVI) were prepared to predict the landslide zonation. Landsat 8 OLI a spatial resolution of with 30 m and ASTER GDEM with a spatial resolution of 30 m (Table 1) were used to generate different thematic maps. ArcGIS, QGIS, ENVI, ERDAS Imagine, PCI GEOMATIC and SPSS software packages were used for data processing.

Methodology
The methodology involves a five-stage analysis: (a) preparation of landslide inventory, (b) selection of geoenvironmental factors, (c) preparation of causative factors, (d) generation of landslide susceptibility map and (e) validation and comparison of the model (Figure 2)

Preparation of landslide inventory
It is necessary to know that the previous landslide-affected area to predict future landslide events (Corominas et al. 2014;Galli et al. 2008;Van Westen 2012). Landslide inventories are the easiest and direct method to map the landslide locations (Guzzetti et al. 2012). The landslide inventory can be verified using the visual interpretation technique of the satellite image, Google Earth and GPS satellite data, historical data and field investigation. Generally, researchers use geomorphologic methods to prepare landslide inventory (Mohammadi et al. 2018;Soeters and Van Westen 1984;Van Westen and Soeters 2006;Youssef et al. 2009). The landslide inventory map at the 1:30,000 scale was prepared using landslide locations based on landsat 8 OLI, Google Earth pro, NASA catalogue and GSI bhukosh and overlaid on the slope map derived from ASTER GDEM ( Figure 3). In this study, a total number of 237 landslides encompassing an area of 943 km 2 were identified; in addition, we randomly collected 237 non-landslide locations, as shown in Figure 3(a). Some landslides digitized from google earth and satellite data are shown (Figure 3(b), Figure 3 (c)). All the landslides were delineated in polygon and converted into point data. Total landslide inventory was then randomly split into two data sets comprising 166 landslides for training (70%) and 71 landslides for testing (30%) using ArcGIS 10.4.

Selection of geo-environmental factors 3.1.2.1 Multicollinearity analysis.
It is necessary to select the proper landslide influencing factors for landslide susceptibility mapping. In the present study, 12 factors were chosen for multicollinearity analysis. Multicollinearity analysis is one of the important principle hypotheses of statistical models. The multicollinearity test enhances the results of the models by selecting the landslide factors as independent factors and landslide inventory points as a dependent factor (Zhang et al. 2018). Multicollinearity is used for various purposes, like piping erosion (Chen et al. 2021b), groundwater potential mapping (Chen et al. 2021b) and flood mapping (Aljuaidi's et al., 2018). In the present study, multicollinearity was tested using the variance inflation factor (VIF) and tolerance (TOL) (Zhang et al. 2018;Lee et al. 2018). The VIF and tolerance value were calculated using equations (1) and (2), where R 2 j is the coefficient of regression of vector j. The coefficient of tolerance less than 0.1 and VIF greater than 5 show the multicollinearity problem (O'Brien 2007). This research tested the multicollinearity of twelve factors using SPSS software.

Preparation of causative factors
For the LSM preparation, there is no universal guideline considered for selecting the factors influencing landslides in this study. It is essential to presume that future landslides may occur at the previous landslide locations (Rasyid et al. 2016;Du et al. 2017). Causative factors selected in this study are topographic, hydrological, geological and anthropogenic factors (Basu and Pal 2020;Nguyen et al. 2019).

Topographical parameters.
Topographical factors like slope, aspect, curvature, elevation and TWI were generated using ASTER GDEM of a spatial resolution of 30 m. An elevation map is one of the most important parameters for landslides; it shows the high concentration of slope instability with increased elevation (Figure 4(a)). Erosion, rainfall/ hydrological conditions and weathering are changed with elevation. The elevation map ranges from 272 m to 4669 m and is classified into 10 classes using ArcGIS 10.4 (Dikshit et al. 2020;Pham et al. 2020;Yilmaz, Topal, and Süzen 2012). The slope is an important factor in landslide mapping. The slope factor and landslide show a positive relationship (Dikshit et al. 2020;Du et al. 2017;Lee and Min 2001). As the slope increases, there are more chances of a landslide because the slope is directly associated with the soil shear stress or other incoherent material increases. The higher slope leads to an increase in shear stress and reduces shear strength, i.e. higher slope areas are more pronounced to landslide than the lower ones. The slope was computed from the ASTER GDEM of 30 m resolution and classified into 10 classes, ranging from 0° to 70° (Figure 4(b)). The aspect was adapted to a compass direction of slope, which varies from 0° to 360°. The slope aspect highly influences weather conditions ( (Figure 4(c)). One of the important topographic attributes is the considered curvature; it shows the morphology of topography. The curvature was calculated from AST5ER GDEM of 30 m resolution and classified into three classes: convex, concave and flat (Hong et al. 2018) (Figure 4(d)). The positive curvature indicates the convex at the grid cell, the negative curvature indicates the concave at the grid cell and value zero indicates the flat surface. The concave curvature has the moisture-holding capacity for a longer period, which reduces the soil's cohesion, whilst the convex curvature is exposed to the expansion and contraction process. More positive and negative curvatures lead to a high drainage concentration, slope instability and slope saturation (Pradhan and Kim 2017; Chen and Li 2020; Kaur et al. 2019).

Hydrological factors.
Distance from the drainage map, also called the drainage buffer map, is a controlling factor for the stability of the slope, soil erosion and the degree of saturation (Akgun et al. 2008;Mishra and Sarkar 2020;Mukherjee et al. 2020). The drainage distance map was prepared using the euclidean distance tool in ArcGIS and classified into ten classes ranging from 0 to 1797.25 m (Figure 4(e)). Rainfall is one of the triggering factors for landslides (Chen et al. 2014;Dikshit and Satyam 2019;Teja et. al. 2019). The increase in precipitation leads to more landslides because the soil becomes saturated, the water level rises and it creates soil erosion, which causes more landslides (Jordanova et al. 2020;Tran et al. 2019). The rainfall map was prepared using CHIRPS satellite data and classified into ten classes ranging from 799.78-1654.35 mm (Figure 4(f)). TWI is the topographic wetness index (TWI). It shows the complexity of the topographic water surface and the flow accumulation tendency in the water basin (Althuwaynee et al. 2014;Moore et al. 1991;Rózycka et al. 2017). The TWI map was prepared based on 30 m ASTER GDEM data using equation 3. TWI ranges from 2.42 to 23.28 (Figure 4(g)),

Geology and distance from lineaments.
Geology plays a vital role in slope instability. Weathered rocks are more vulnerable than hard rocks. The structure of geological features affects the sliding of rocks (Dai et al. 2002;Fell et al. 2008; L. M. Highland and Bobrowsky 2008). In this study, the geology map was digitized on the scale of 1:30,000 collected from the geological survey of India and classified into four groups: Kanchenjunga gneiss, Chungthang schists and gneiss, Darjeeling gneiss and Daling (Figure 4(h)). These groups belong to Precambrian age. Chungthang schists and gneiss include quartz-biotite schist, Al-silicate rocks and graphite schists. Kanchenjunga rocks are hard, dense, well joined and intruded at places by pegmatite and tourmaline granite (Anbalagan 1996). Darjeeling Gneiss consists of genetic group rock from Precambrian to the late Tertiary times. The Darjeeling group accommodated isolated lenticular and twisted, fine-grained, concentric calc-silicate eye rocks. The Daling represent schistose; they are composed of phyllites, sericite and chlorite-phyllites, slates, flagstone, quartzite, greywackes, conglomerates, silt-stone and schists (Kaur et al. 2019). Two significant thrusts in the East Sikkim are the Chungthang Thrust and the Main Central Thrust. Different rock units experience variable tectonic activities and exhibit high weathering near the drainage channel and the thrust/fault (Anbalagan et al. 2015;. Lineament is a surficial expression of subsurface geological features such as fault, joint and fracture (Pal et al. 2006(Pal et al. , 2007Narayan et al. 2017). Lineament represents the weaker zone, and almost all  Jaafari et al. 2015). Most of the landslide occurs near the road. The roads were digitized from the GSI Bhukosh. The distance from the road was the estimated and divided into ten buffer classes ranging from 0-12,668 m using the Euclidean distance tool in ArcGIS (Figure 4(j)).

Normalized Difference Vegetation Index (NDVI) and Land Cover Land Use (LULC). NDVI was
used to identify the vegetation on the surface and was calculated using the formula NDVI = (IR-R)/(IR+R), where IR represents the infrared band and R represents the red band (Ba et al. 2018;Behling et al. 2014;Mishra and Sarkar 2020). The NDVI map was prepared using Landsat 8 OLI (Figure 4(k)). The NDVI value varies from −1 to +1, where the negative value indicates the area of non-vegetation and the positive value shows the vegetation area. The LULC map was prepared from Landsat 8 OLI of the 2019 satellite image based on supervised classification techniques in ERDAS Imagine software, which was also verified with BHUVAN and Google Earth data. All the classes were digitized with the help of visual interpretation and NDVI of Landsat 8. The LULC map was classified into six classes, i.e. settlement, agriculture land, barren land, medium forest, thick forest and spare forest (Figure 4(l)). Classification accuracy indicated that the kappa coefficient was 80.5%, which predicted good accuracy. The classified map showed the settlement of 14.021%, an agriculture land of 3.003%, a medium forest of 28.45%, a thick forest of 27.43%, a spare forest of 12.47% and a barren land of 14.81%

Frequency ratio (FR) method
The frequency ratio (FR) method was proposed by Lee and Talib (2005). FR is the ratio of percentage landslide in the class to the percentage of the class area of the entire map (Lee and Pradhan 2007). A major reason to utilize this method is the basic simplicity of the methods and quick process of a large amount of data in a typical GIS environment. The frequency ratio of individual pixels where S pix (Li) is the pixel count having landslide in each pixel in every class, S pix (Ni) is the total pixel count in the area having class (i) = ∑S (Li), i.e. total pixel count in an individual class, and ∑S pix (Ni) is the total pixel count in the area. This method is performed using GIS software. Then, the calculated FR model of individual factor value was summed in ArcGIS to calculate the LSM, FR indicates a higher correlation for a value greater than 1, whilst a value less than 1 indicates a lower correlation in the area.

Logistic regression (LR)
Logistic regression (LR) is a multivariate regression model that shows a relationship between the dependent and independent variables (Goetz et al. 2015;Martha et al. 2017). The advantage of using LR over simple regression is that by adding a congruous function to the usual linear regression models, variables may be categorical or continuous or a combination of both categorical and continuous. In the multivariate regression, the dependent variable factors must be binary, whilst the bivariate statistical model and discriminate analysis variables must have a normal distribution. The algorithm applies in the logistic regression model maximum likelihood estimation after converting the dependent variable into 0 or 1 or logit variable (Cama et al. 2016;Pourghasemi et al. 2018). In landslide susceptibility mapping, the dependent variables are binary numbers representing the presence (1) or absence of landslide (0). In contrast, independent variables are 12 Geoenvironmental factors that influence landslides (Table 1). With the help of dependent and independent variables, LR finds the best fitting line model to show the relationship between the presence and absence of landslide Pradhan 2010). The relationship between dependent and independent variables is non-linear. The logistic regression can be expressed according to equation 6 (Lee and Min 2001) Z ¼ β 0 þ β 1 X 1 þ β 2 X 2 þ ::β n X n ; where p is the probability of the event occurring varying from 0 to 1 and its S-shaped curve, z is the linear combination, β 0 is the intercept of the model, β 1 to β n are the slope coefficient of the logistic regression and the X 1 to X n are the independent variables. A positive sign indicates that the variables positively correlate with changes, and a negative sign implies the opposite effect.

Random forest
The Random Forest (RF) is a supervised machine learning classification method developed by Leo Breiman (2001), which has been considered to have more potential for classification and regression learning. Random forest is a tree-based ensemble technique that uses the bootstrapping or bagging algorithm. The advantage of the RF model is to minimize the overfitting problem by building multiple decision trees, and from the random subset, they split the node on the best split ). Each training sample was built on bootstrap aggregating techniques randomly selected from the training data observation, and then final classification results were voted from all the decision trees.

Analytical hierarchy process with FR (AHP-FR)
The AHP method developed and modified by Saaty (1980Saaty ( ), (1987Saaty ( ), (1994 provides an easy and flexible understanding to solve a complicated problem. AHP is a multicriteria decision method that allows solving the complicated problem in terms of an eigenvalue through the pairwise comparison techniques. The calculation of eigenvalues is part of the AHP model. It is necessary to break the complex problem into the component factors, arrange them in hierarchic order (Yalcin et al. 2011) and then assign the weight to each factor according to the importance. The final stage adapts a rating to overall priorities factors (Saaty, 1994). The AHP requires a scale ranging from 1 to 9 to compute the quantitative and qualitative performance (Table 2), where 1 indicates the lowest importance and 9 expresses the highest priority to the final solution (Saaty, 1994;Saaty 1987). For integrating with the frequency ratio (FR), first, the relative frequency of each class of factor was calculated using equation 8, where RF is the relative frequency and FR ij is the frequency ratio The next stage is to calculate the predictive rate (PR) for each factor, where PR is the the prediction rate, RF imax is the maximum relative frequency of factor i, RF imin is the minimum relative frequency of factor i, ðRF max À RF min Þmin is the minimum difference between max and min relative frequencies of all the factors.
Then, the calculated PR weights were used to estimate the AHP using a pairwise matrix.
AHP requires the creation of a pairwise comparison matrix that is subjective, and results are highly dependent on expert judgments (Ergu et al. 2011). In AHP, the consistency ratio indicated the probability of the randomly generated judgment matrix estimated using equation 10, where RI is the average of the resulting consistency index and CI is the consistency index shown in equation 11, where λ max is the largest principal eigenvalue of the matrix. The model with a CR value greater than 0.1 requires a revision of the judgment, but the CR below the 0.1 model is accepted.

Model validation and accuracy
There are different statistical indexes such as accuracy (ACC), specificity, sensitivity, positive predictive rate (PPR), negative predictive rate (NPR), kappa index, root mean square error (RMSE), success rate and prediction rate curve and seed cell area index (SCAI). Many authors used the ROC for the validation of the model. The ROC is the most common and popular statistical method that has been applied in landslide susceptibility mapping (Hongb et al. 2016;Lei et al. 2021). ROC is important for comparing different landslide susceptibility models, by computing their success rate and prediction rate curves. Performance and validation of landslide mapping were performed under the specificity, sensitivity, PPR, NPR, ACC, RMSE, kappa, SCAI, success rate and prediction rate curve (Chung and Fabbri 2003;Süzen and Doyuran 2004;Yilmaz et al. 2012). Sensitivity is the proportion of landslide pixels that are correctly classified as landslide occurrence (Equation 12). Specificity is the proportion of non-landslide pixels that are correctly classified as non-landslide occurrences (Equation 13). Accuracy is the ratio of landslide pixel and nonlandslide pixels that are correctly classified (Equation 14). The positive prediction rate is the probability of pixels classified as landslide (Equation 15), and the negative prediction rate is the probability of pixels classified as non-landslide (Equation 16) Bui et al. 2016;Chen et al. 2018b). They could be calculated using the following equations: where TN (true negative) and TP (true positive) are the numbers of pixels correctly classified, whilst FP (false positive) and FN (false negative) are the numbers of pixels incorrectly.

ROC curves
The success rate curve is 'goodness of model fit', which was acquired by comparing a landslide-training data set with a landslide susceptibility map (HembRam et al.  2020; Xu et al. 2013). The area under the curve (AUC) of prediction rate curves was counted based on each landslide susceptibility map and testing landslide data set. The landslide susceptibility map was reclassified into 100 equal classes and sorted into descending order to calculate the success rate curve and prediction rate curve. The tabulated area tool in ArcGIS was used to calculate the landslide pixels in each susceptible class using the training data set to calculate the success rate curve and testing data for the prediction rate curve. The cumulative percentage was calculated, and then the graph was  the high accuracy of the model, whilst a high SCAI value indicates the poor accuracy of the model. The value (Süzen and Doyuran 2004) give the SCAI formula and was calculated as follows:

Multicollinearity analysis of selected independent variables
The multicollinearity amongst independent factors was identified using tolerance (TOL) and variance factors (VIF). The results show that each independent variable has VIF<5 and TOL >0.1. The highest multicollinearity VIF value is 4.825 for elevation, which is less than 5. Therefore, all the independent variables need to be selected for landslide susceptibility mapping. The TOL and VIF of each independent parameter are shown in Table 3. The multicollinearity analysis results indicates that there is no collinearity problem between the independent variable set.

Landslide susceptibility using FR models
The landslide susceptibility map has been prepared by integrating the entire thematic layer with the weights using the FR method. By analysing the relationship of 12 geo-environmental factors with the landslide training occurrence, the FR ratio was calculated (Table 4). The pixel number of every class of causative factors was calculated with the help of a reclassification tool in ArcGIS. It can be seen from Table 4 that the relationship of landslide occurrence with the elevation of higher FR weight (2.83) occurred in the range of 2638 m -3181 m. It was clearly shown that at a higher elevation, the FR value decreases. At an elevation of 4022 m-4669 m, landslides were rarely identified and exhibited a significantly less relationship at the higher elevation. The slope and landslide exhibit a positive relationship; a higher slope promotes more chances of a landslide. A slope greater than 30° has an FR value greater than 1, indicating a higher chance of landslide at the higher slope. At the highest slope, the FR value is 5.36, indicating a greater chance of landslide, whilst at the lowest, the FR value is 0.80, indicating a lesser chance of landslide occurrence. From the aspect map, the higher FR values of 3.52 and 2.42 were estimated in the south and southeast directions, respectively, which indicates the number of landslide occurrences of those directions. The south-facing slope witnessed higher solar heating, rainfall and weathering than the north-facing. It shows that the south and southeast-facing slopes are more unstable. In curvature, concave a higher FR value of 1.36 was found compared to convex because the concave surface retains moisture due to heavy rainfall for a long time. Due to weathering, the soil becomes weak, causing a higher FR value showing a greater landslide probability. The TWI shows that the strength of soils decreases with the increase in pore water pressure. In the study, TWI class 15.36-23.28 has the highest FR value of 3.92. The rainfall class ranging from 1148.81-1181.82 mm exhibits a higher FR value of 1.76, whilst a lower amount of rainfall exhibits a lower FR value of 0.41. The distance from drainage is directly related to the rainfall pattern in the area; the FR value decreases as we go to higher drainage density; however, the probability of landslide occurrence is relatively high near the drainage. A higher FR value of 2.47 was found near the stream for class 0-133.913 m, and the area away from the stream indicated a lower FR value. Considering the distance from the lineament, the highest probability of landslide occurred near the lineament having a FR value of 1.45 found between 0 and 258.16 m, whilst as we go higher, the FR value decreases.  In the case of geology, chungthang schists and gneiss group rock was more susceptible to landslide, having an FR value of 1.641, followed by Darjeeling gneiss having an FR value of 0.955. Chungthang schists and gneiss consisting of quartz, biotite schist and slate are more susceptible to landslide due to the presence of MCT in that area. Considering the LULC, a higher FR value was found in settlement (FR = 3.99) and agriculture (FR = 2.56), indicating more landslide vulnerability. The higher landslide vulnerability of the settlement area means that human activity affects the landslide occurrence. The highest FR value of 1.487 was found in the NDVI class 0.17-0.21, indicating that lower vegetation cover area is more susceptible to a landslide. A higher landslide probability was found near the road with a higher FR value of 2.04.
The final landslide susceptibility map was organized by summing all the causative factors (Table 4) using the following equation in the raster calculator function of the spatial analyst tool in ArcGIS 10.4. The generated landslide susceptible map (LSM) was reclassified into five groups using the natural break method. LSM represents the probability of landslide occurrence in the future. The greater the value, the higher the landslide vulnerability, The LSM map of the East Sikkim area was divided into five categories, i.e. very low, low, moderate, high and very high, as shown in Figure 5 based on the natural break method, which covers 20.74%, 29.80%, 25.56%, 16.75% and 7.13% area, respectively (Figure 9(d)). The prepared landslide susceptible map shows that 7.13 % of landslide comes under the very high vulnerable zone, whilst 20.74 % exhibit a low vulnerability zone (Table 5).

Landslide susceptibility using logistic regression (LR)
The landslide susceptibility map was prepared using logistic regression (LR) in python, based on the dependent and independent factors. Initially, all the factors were converted into ASCII format. LR coefficients of curvature, slope, TWI, elevation, geology and LULC show positive correlation, whilst aspect, distance from drainage, distance from lineament and rainfall have negative correlations (Table 6). All the calculated values show that all factors were equally important for the occurrence of landslides. Then, final LSM was prepared by using equation (20), which predicts landslide occurrence possibilities, LR ¼ ðslopex0:069Þnðaspectx0:006Þ þ ðcurvaturex0:103Þ þ ðTWIx0:097Þ þ ðelevationx:038Þnðdistance from drainagex:001Þ À ðrainfallx0:005Þ þ ðdistance from roadx0Þnðdistance from lineamentx:002Þþ ðgeologyx0:952Þ þ ðLULCx:103Þ nðNDVIx1:037Þ þ 8:586: The prepared LSM value lies between 0 and 1, has the pixel size of 30 m x 30 m and reclassified into five classes as very low, low, moderate, high and very high shown in Figure 6, which cover the area of 10.90%, 22.44%, 28.21%, 26.44% and 11.99 % (Figure 9(b)), respectively, as shown in Table 5.

Landslide susceptibility using AHP-FR
The AHP analysis indicated the CR value of 0.05, exhibiting a good level of consistency in the pairwise matrix, and it is suitable for assigning weights to each factor. Using the AHP method, weights of each factor were calculated based on the FR method (Table 7). According to Table 7, slope was the most influencing factor on a landslide with the AHP-FR weight of 0.14 and then the LULC with the AHP-FR weight of 0.12. However, the curvature, NDVI, and geology exhibited lower AHP weights of 0.01, 0.03 and 0.03, respectively, whilst aspect, distance from drainage, distance from the road, rainfall, distance from lineaments and TWI exhibited higher AHP weights of 0.10, 0.07, 0.06, 0.05, 0.04 and 0.1, respectively. An integrated AHP-FR-based LSM map was produced by using the following equation: The LSM value was then reclassified into five classes as very low, low, moderate, high and very high classes (Figure 7), covering 30.05%, 23.37%, 18.86%, 15.73% and 11.97% of the area (Figure 9(a)), respectively (Table 5).

Landslide susceptibility mapping using random forest
The landslide susceptibility map (LSM) was prepared using the Random Forest (RF) model and classified into five classes (very high, high, moderate, low and very low) using the natural break method (Figure 8). The very high landslide susceptibility zone occupied about 12.69 % of the East Sikkim district, whilst the high, moderate, low and very low zones occupied 21.63%, 27.75%, 12.98% and 24.95%, respectively (Figure 9(c)). Approximately 73.73% of landslide pixels were observed in the very high zone. Lithologically, this zone is highly vulnerable to landslide occurrences.

Validation and comparison of landslide susceptibility models
Models were validated using different statistical methods such as sensitivity, ACC, NPR, PPR, specificity, RMSE, Kappa, success rate and prediction rate curve for both the training and testing data sets of AUC curves. The RF model shows the best performance for the training data set, followed by the LR model. In the case of ACC, RF had the highest accuracy of 81.02 %, followed by the LR model, 73.19%. In the case of sensitivity and specificity, the RF model had the 82.53% and 79.51%, followed by the LR model 72.51% and 73.91%. In the case of PPR and NPR, the RF model had 80.11% and 81.98% followed by the LR model, 74.69% and 71.68%. In addition, the kappa value of the RF model is 0.61 and 0.46 for the LR model. The RMSE value of the RF model is 0.37 and 0.44 for the LR model (Table 8a).
For the validation data set, in the case of ACC, RF had the highest accuracy of 76.76% followed by the LR model, 72.53%. In the case of sensitivity and specificity, the RF model had 77.46% and 76.05% followed by the LR model, 77.38% and 67.60%. In the case of PPR and NPR, the RF model had 76.38 % and 77.14% followed by the LR model, 70.51% and 75.00%. In addition, the kappa value of the RF model is 0.53 and 0.45 for the LR model. The RMSE value of the RF model is 0.38 and 0.42 for the LR model (Table 8b).
The success rate curve and prediction rate curve of landslide susceptibility maps obtained from AHP-FR, RF, LR and FR models are shown in Figure 10. The prediction rate curve shows how well the landslide models and landslide factors are performed. The AUC of the success rate curve exhibited 0.88 for the RF model, 0.85 for AHP-FR, 0.78 for LR and 0.79 for FR, as shown in Figure 10(a). The prediction rate curve exhibited 0.86 for the RF model, 0.81 for AHP-FR, 0.79 for LR and 0.78 for FR, as shown in Figure 10(b). An AUC equal to 0.5 indicates the least accuracy, whilst an AUC equal to 1 indicates the best prediction accuracy. The RF method indicated the best performance on both training and testing data sets having the highest accuracy of 0.88 for the success rate curve and 0.86 for the prediction rate curve, followed by the performance of AHP-FR for both success and prediction rate curves. This indicated that machine learning was the best result and also showed that the integration of heuristic models with statistical modes indicating a better result than the other statistical models.The LR and FR methods have relatively slightly lower AUC.
Overall results revealed that all models have good prediction capabilities. The validation of LSM classes and landslide pixel of testing were computed using the SCAI statistical methods for four models. From Table 9, it is observed that the value of SCAI decreased from very low classes to very high classes. The SCAI values of very high susceptibility classes are 0.14, 0.17, 0.18 and 0.19 for the RF, AHP-FR, LR and FR models, respectively. The result of the SCAI method was perfectly matched with the rules, so the accuracy and prediction of all the models were very satisfactory.

Discussion
Landslide is one of the most hazardous natural disasters in mountainous areas that can cause a huge loss of lives and property, leading to massive socioeconomic and environmental effects. The highways and settlements witnessed high vulnerability owing to the dynamic nature and unconsolidated bedrocks of the East Sikkim. This problem mainly occurs in the monsoon season because of heavy rain causing increased pore water pressure. In the last two decades, the landslide has increased due to changes in land use land cover and increase in unplanned settlements, deforestation and infrastructure development.
Previous studies investigated the landslide locations covering a small part of East Sikkim using the knowledge-based and statistical method. The present study considered the whole East Sikkim region using knowledge-based (AHP-FR), statistical (FR and LR) and machine learning (RF) approaches. All models are popular and widely used methods for landslide mapping and have their specific advantage and disadvantage. The heuristic model (AHP) involves knowledge-based analysis in which weights are assigned based on expert opinion and field investigation to evaluate the importance of landslide causative factors. But sometimes, knowledgebased (subjective) methods fail to provide reliable results due to different opinions of the experts that lie to determine the uncertainty produced by the changes in the weight and uncertainty of the factors that would affect the final landslide susceptibility maps. Both the statistical methods (FR and LR) are based on objective analysis in which the model computes weights and rating to each factor class based on the logic tool. But statistical models require more data to generate reliable results, so machine learning methods are widely used to overcome the performance of the statistical method.
Generally, heuristic methods are time-consuming and expert knowledge does not provide an accurate result. Therefore, uncertainty associated with the weights of factors for each class is minimized by integrating the heuristic method with the statistical method. Then, we finally compare the performance of the objective (statistical) and subjective (knowledge-based) methods with that of the machine learning method.
(1) The LSM prepared using the AHP-FR, LR, RF and FR models was classified into five classes, very low, low, moderate, high and very high using natural break classification methods in ArcGIS. In East Sikkim, very high classes of the RF, AHP-FR, LR and FR models cover the area of 12.69%, 11.97%, 11.99% and 7.13%. After comparing the performance of different models using sensitivity, specificity, ACC, NPR, RMSE and kappa index, it is found that the RF model shows the best result for both training and testing data. The ROC curve of success rate and prediction rate curve for the RF model have higher prediction performance (88 %) followed by AHP-FR (85%) than logistic regression (79%) and frequency ratio models (78%). In the case of the prediction rate curve using testing variables, 86% for RF, 81% for AHP-FR, 79% for LR and 78% for FR were observed. However, the study shows that machine learning performed well, whilst the performance of the heuristic model can be improved by integrating with a statistical model.

Conclusion
In the present study, twelve landslide susceptibility thematic maps were prepared in ArcGIS software using slope, aspect, elevation, curvature, distance from lineament, distance from drainage, distance from road, LULC, NDVI, TWI, geology and rainfall. These thematic maps were generated using the Landsat 8 OLI and ASTER GDEM data. The obtained results from RF, AHP-FR, LR and FR susceptibility maps indicated that the RF and AHP-FR model had a higher success rate and prediction rate curve. In summary, we can conclude that the integration of heuristic classifiers with the statistical model improves the results of the heuristic method (AHP). This study has some limitations; we could not include soilrelated data in this study. If the surrounding area is recognized as relatively highly prone, a more detailed geotechnical study should be carried out to assess the risk further before construction or development begins for the accurate susceptibility map. More machine learning concepts can also be used for better results in landslide susceptibility assessment. The hazard susceptibility categories interpreted both spatially and quantitatively shall be important for further development of infrastructures like roads and settlements in the area. Required planning level precautions and mitigation efforts can be    suitably supplemented to minimize the incidences of landslide hazards and environmental degradation with these Landslide Hazard Zonation (LHZ) maps produced in the medium scale for the area. Continuous updates and an increase in the accuracy of these maps could help prevent huge damage and mitigate loss.