Evaluating the Performance of Individual and Novel Ensemble of Machine Learning and Statistical Models for Landslide Susceptibility Assessment at Rudraprayag District of Garhwal Himalaya

: Landslides are known as the world’s most dangerous threat in mountainous regions and pose a critical obstacle for both economic and infrastructural progress. It is, therefore, quite relevant to discuss the pattern of spatial incidence of this phenomenon. The current research manifests a set of individual and ensemble of machine learning and probabilistic approaches like an artificial neural network (ANN), support vector machine (SVM), random forest (RF), logistic regression (LR), and their ensembles such as ANN ‐ RF, ANN ‐ SVM, SVM ‐ RF, SVM ‐ LR, LR ‐ RF, LR ‐ ANN, ANN ‐ LR ‐ RF, ANN ‐ RF ‐ SVM, ANN ‐ SVM ‐ LR, RF ‐ SVM ‐ LR, and ANN ‐ RF ‐ SVM ‐ LR for mapping landslide susceptibility in Rudraprayag district of Garhwal Himalaya, India. A landslide inventory map along with sixteen landslide conditioning factors (LCFs) was used. Randomly partitioned sets of 70%:30% were used to ascertain the goodness of fit and predictive ability of the models. The contribution of LCFs was analyzed using the RF model. The altitude and drainage density were found to be the responsible factors in causing the landslide in the study area according to the RF model. The robustness of models was assessed through three threshold dependent measures, i.e., receiver operating characteristic (ROC), precision and accuracy, and two threshold independent measures, i.e., mean ‐ absolute ‐ error (MAE) and root ‐ mean ‐ square ‐ error (RMSE). Finally, using the compound factor (CF) method, the models were prioritized based on the results of the validation methods to choose best model. Results show that ANN ‐ RF ‐ LR indicated a realistic finding, concentrating only on 17.74% of the study area as highly susceptible to landslide. The ANN ‐ RF ‐ LR ensemble demonstrated the highest goodness of fit and predictive capacity with respective values of 87.83% (area under the success rate curve) and 93.98% (area under prediction rate curve), and the highest robustness correspondingly. These attempts will play a significant role in ensemble modeling, in building reliable and comprehensive models. The proposed ANN ‐ RF ‐ LR ensemble model may be used in the other geographic areas having similar geo ‐ environmental conditions. It may also be used in other types of geo ‐ hazard modeling. of spatial relationship between causative factors and previous landslides. In this work, we used landslide causative factors using remote sensing data and applied several machine learning ROC curve, efficiency, accuracy, MAE, and RMSE were employed to assess and compare these LSS models. Finally, using the CF the best model was selected. The findings of the validation techniques show that the LSS models prepared by ANN, FR, SVM, LR, and their ensembles performed well in this analysis. The ensemble of ANN, RF, SVM, and LR models played an improving role as they improved the values of AUC. However, with respect to the results of robustness, the ANN ‐ RF ‐ LR ensemble had the highest reliability in terms of results where the highest agreement was found in the efficiency, accuracy, AUC of success rate curve, and predictive rate curve values, MAE and RMSE. The present research shows that ensemble models are effective tools that could be used in landslide prediction and landslide management for efficient decision ‐ making by the authorities. In the present analysis, these models were utilized in the Rudraprayag district. However, for wider use, the output of such models may be used for areas of similar geo ‐ environmental conditions. The key drawback of this analysis is that we used a limited set of LCFs and a fixed training testing samples dataset ratio (70:30) for the all the models. The performance of the models will have to be further tested through a variety of training and testing sets of sample datasets as well as LCF sets.


Introduction
The landmass movement is the physical process in downhill geomorphic regions, which causes rolls of rubble, regoliths, and large masses of soil in down slope direction under the impact of gravity. It is a predominantly geological event which takes place when the force of the material exceeds the resistance of the shear strength of the soil [1]. Landslides are ranked as the seventh most destructive geohazard among the various natural hazards prevalent around the world in terms of the magnitude of the loss of lives and properties [2]. The occurrence of landslides depends upon the physiographic, hydrological, geological, and geomorphic setups of the study area. Sometimes, the movement of a large landmass in the hilly regions causes a change in topography, drainage pattern, and courses of rivers. The impact of landslides can be attributed to a wide range of phenomenon such as deforestation, habitat loss of human and wild lives, soil loss, disruption of roads and constructions, and flooding in specific cases [3,4]. In the mountainous landscape, landslides occur frequently, triggered by a heavy rainstorm along with interactions among the various environmental causative factors. According to Raman and Punia [5], around 0.49 million sq.km land of India is under the threat of landslide hazard, which is about 15% of the land area of the country. In Indian sub-continent, landmass movement and failure of slope are mostly activated by a heavy precipitation during monsoon, however in the Himalayan belt, the occurrence of landfalls are pretty normal due to seismic activities viz. the Uttarakhand earthquake in 2017, the Chamoli earthquake in 1999, and the Sikkim earthquake in 2011 [6][7][8][9]. Wadhawan [10] mentioned that the mountainous region of Konkan and the Nilgiri ranges in the south and south-western India and the Western Ghats are severely vulnerable to landslides and mass movements. Other than seismic tremors and rainstorm, human interventions and unplanned encroachment, i.e., construction of roads, civil structures etc., disturb the resistance of landmass and accelerate the rate of landslide occurrence.
Therefore, it is important to know the spatial relationship between the geo-environmental factors and landslide events for quantitative evaluation of landslide susceptibility (LSS). The reliability of analysis and modelling is determined by the quality and availability of input data and methodology employed. Assembling landslide data in mountainous regions is a very challenging task, mostly owing to poor accessibility. Obtaining data through field inspection techniques can be time intensive and costly, and application of remotely sensed data in this regard is fruitful. The literatures suggested that landslide susceptibility zonation has been performed previously using various bivariate and multivariate statistical techniques, such as frequency ratio [11,12], evidential belief function [13,14], statistical index model [15], logistic regression [12,16], weight-of-evidence [15,17], etc. These studies have shown a good accuracy of these techniques in predicting LSS spatially. However, with the advancement of the machine learning-based modelling, LSS has come out as a greater efficient and effective approach in recent decades. There are various studies that used machine learning and artificial intelligence models in performing landslide susceptibility (LSS) assessment such as random forest [18,19], boosted regression tree [12], naïve bayes [20], decision tree [21], neuro-fuzzy [14], artificial neural network [22], support vector machine [20], and many others.
However, one of the major weaknesses of statistical models lies in the fact that some assumptions should be defined before conducting any analysis, and therefore a spatial relationship between causative factors is widely neglected [23]. In contrast, while dealing with the machine learning models, the advantages are they do not need statistical assumptions before analysis and they can deal with data of various types, i.e., categorical and continuous [13]. In this regard, the ensemble models can resolve nonlinear difficulties and complexities [24]. In recent decades, numerous studies were conducted using the ensemble model approach which has shown greater effectiveness in the spatial assessment of various environmental hazards viz. gully erosion susceptibility assessment [25,26], groundwater contamination risk assessment [27], groundwater potential mapping [28], flood risk assessment [29], and landslide vulnerability assessment [30,31].
The main purpose of the predictive models is to recognize the LSS extents depending upon the nature of spatial relationship between causative factors and previous landslides. In this work, we used landslide causative factors using remote sensing data and applied several machine learning models to produce landslide susceptibility maps (LSMs) of Rudraprayag district in Uttarakhand. Subsequently, LSMs were ensemble to attain better accuracy that can be achieved from a single model. The ensemble of the models in the present analysis was implemented in four-stages, whereby in the first stage, individual models were calibrated and then a two-model combination approach was performed, i.e., artificial neural network-random forest (ANN-RF), ANN-support vector machine (ANN-SVM), SVM-RF, SVM-logistic regression (SVM-LR), LR-RF, and LR-ANN. In the second stage, LSS model calibration was performed compiling three model ensemble approaches, i.e., ANN-LR-RF, ANN-RF-SVM, ANN-SVM-LR, RF-SVM-LR. Finally, all four models were combined.

Study Area
The Rudraprayag district is located in a highly elevated (2936 feet or 895 m) part of the Garhwal Himalaya in the state of Uttarakhand, India ( Figure 1). This downhill district covers an area of 2439 sq. km and geographically extended between 78°49ʹ E and 79°21ʹ2ʹʹ E and from 30°12ʹ N to 30°48ʹ N. According to Hindu mythology, the name of the district comes from the terminology "Rudra" which means "god of wind or storm" in Hindu religion. Because of the unique physiographic set-up of this Himalayan region, the Rudraprayag district is frequently affected by landslides in the past few decades causing deaths and economic losses [32]. The hilly parts of this district are intensely deepened by the major river Alaknanada. Relative relief per unit area is very high (500-3000 m) and the district is categorised under steep-sloped topography. Soil formation is very weak and the layer depth ranges between 0 and 10 m. The geological formation of this area is characterised by Guinness and Dalling sequence, Siwalik rock, Gondwana Park rock and Pleistocene rises. The area lies within the sub-tropical humid region where the highest temperature (25 to 37 °C) prevail in May to July and lowest temperature (2 to 16 °C) recorded between November to February. The average annual precipitation is about 405.17 cm of which 75% to 80% occurs in monsoon (between June and September) when numerous landslides take place (Indian Meteorological Department [33]. Because of its location in a seismic-tectonic disturbance zone and due to anthropogenic interventions (often found along the National Highway 109 and 58 of India), landslides have become a frequent phenomenon in this region. Such topographic, climatic, and anthropogenic frameworks are intensifying the need for the identification of potential areas of landslides in the Rudraprayag district.

Materials Used
Topographical map of the study area (sheet no. 78B/5) in 1: 50,000 scale was collected from the Survey of India (SOI), Kolkata which was visually compiled with Google Earth (GE) images from 2017 to 2019 for updating the river map, settlement, and road map. Land use/land covers categories were classified using Sentinel-2 satellite imagery of 10 × 10 m resolution downloaded from the US Geological Survey Earthexplorer [34] which was crosschecked using Cohen's Kappa coefficient through field truth, GE imagery. The Geological Survey of India provided a detailed lithological map in 1:250,000 scale, which was resampled to match to scale with other thematic layers. To map the seismically disturbed zones, an earthquake zoning map was downloaded and extracted from National Center for Seismology, under the Ministry of Earth Sciences, India in 30 m × 30 m resolution. The lineament was extracted in 1:50,000 scale from the web platform of ISRO (Indian Space Research Organization) that is Bhuvan [35]. Slope, aspect, topographical wetness index (TWI), and stream power index (SPI), were extracted from advanced land observation satellite (ALOS) phased array type L-band synthetic aperture radar (PALSAR) derived digital elevation model (DEM) of 12.5 × 12.5 m resolution downloaded from the Alaska Satellite Facility. Spatial variation of soil depth and texture at 1:50,000 scale was gathered from NBSSLUP Regional Centre, Kolkata and classified based on the texture classification method of U.S. Department of Agriculture (USDA). The rainfall data was collected from the Indian Meteorological Department (Table 1). However, the spatial resolution of these factors are not same. For producing the landslide susceptibility map, the resolution of the ALOS PALSAR DEM (12.5 m × 12.5 m) was considered as the base resolution and all the causative factors were resembled into 12.5 × 12.5 m resolution.

Methodology
This study involved several steps and processes for creating landslide susceptibility (LSS) maps including ( Figure 2): (i) A total of 223 landslide locations were identified using the high-resolution Google earth images and afterward these locations were verified through field investigation with a global positioning system (GPS) which was conducted during April 2018 and September 2019. The same number of non landslide points as landslide locations were taken randomly for training the models. The 16 environmental factors were considered for modeling (Table 1). (ii) Relief-F technique was used to judge the effectiveness of the landslide conditioning factors (LCFs) for LSS mapping. (iii) LSS maps first were prepared using ANN, SVM, LR, and RF models. The ensemble models were prepared combining the two, three and four models subsequently. (iv) The contribution of the LCFs was assessed using the random forest (RF) model, (v) The LSS model's performances were evaluated through the area under receiver operating characteristic curve (AUCROC), precision, accuracy, mean-absolute-error (MAE), and rootmean-square-error (RMSE). (vi) Finally, compound factor (CF) method was used to choose the best model.

Generation of Landslide Inventory (GLI)
Generation of Landslide Inventory (GLI) is the set of previous and present locations of landslides and is the prerequisite for modeling the landslide susceptibility (LSS). LSS modeling can be assessed through understanding the connections between GLI and landslide causative factors [36]. GLI data set can be prepared through compiling various data sources considering field investigation, past landslide records and evaluation of the Google Earth imageries [37]. In this research, the GLI was prepared by compiling a total of 223 landslide locations recorded in between April 2018 and September 2019 through several field campaigns ( Figure 3) using GPS and secondary historical record (record of National Disaster Management Authority, Uttarakhand). Among these inventory, over 40% landslides are translational slides. According to the Uttarakhand National Disaster Management Authority, huge landslides have occurred because of cloudburst rainfall that occurred in July 2019. For training the models, we have selected same amount of non-landslide points randomly over the study area. These historical landslides were then spatially mapped in a 12.5 × 12.5 m resolution, dividing randomly into two sub-set viz. training sub-set, which is comprised of 70% (156 nos.) of GLI and validation or testing sub-set, which includes remaining 30% (67 nos.). The separation of the whole GLI was performed using the extension titled "Geostatistical Analyst" in ArcGIS version of 10.3.1. Training sub-set was used for simulation of the models to generate LSMs and validation sub-set was applied to assess the models' performances ( Figure 1). Landslides mainly occurred in this study area due to rainfall and earthquakes. Rate of landslide has increased with the anthropogenic activities as well. Among the recorded landslides, smallest and largest was 21.73 sq. m and 633.78 sq. m, respectively, while the mean extent was 272.41 sq. m.

Relief-F Method
Research relief-F is an effective method for the selection of the relevant parameters. Kononenko (1994) [38] suggested various relief models. Investigational outcomes have displayed that Relief-F attains greater performance than the original Relief [39]. It synthesizes the factor weights following their relevance with the goal. Suppose 'P' is a randomly selected sample. Then, 'RWi' (Relief-F weight) of i-th item can be assessed following Equation (1). One random point can estimate its two nearest neighbors point (one is its same class; other is from diverse class). ← here, NH = nearest hit and NM = nearest miss (1)

Preparation of the Landslide Causative Factors (LCFs)
Selection of the causative factors in an analysis is very challenging as it has no determined norms. Similar hydrological, topographical, and environmental variables are not appropriate for modeling and may be varied according to the geographic setting of the regions. In this research, the LCFs were chosen based on the existing literature [10,40,41] and applying the one factor selection method, i.e., Relief-F. For this study, sixteen LCFs were selected for landslide susceptibility mapping. Among the selected factors, rainfall and earthquake zone are landslide triggering factors and the rest are landslide causative factors. The details of the LCFs are given in Table 1.
Here, a to i indicates the cell value of 3 × 3 window. [44] Curvatu re where Z1-Z9 are altitude values in 3 × 3 cellular networks and S  denotes the cell size.
[45] SPI tan( ) SPI As slope   where 'As' is the specific catchment area in meters. [ where, NIR is near infrared band and IR is the infrared band. [51] LU/LC Supervised classification (Maximum likelihood) [10] Topographic LCFs were extracted from a PALSAR DEM with the resolution of 12.5 × 12.5 m using 'surface' tool in ArcGIS 10.3.1 software. These LCFs are altitude, slope gradient, slope aspect, curvature and stream power index. Higher altitude areas, especially downhill segments of hilly regions, are favourable for landslides (LS) and the study area has a range of 500 ft. to 2936 ft. elevation ( Figure 4a). Slope gradient (Figure 4b) is the result of a combination of factors, i.e., drainage intensity, vegetation, run-off, and geologic. Thus, it was selected as one of the key LCFs in studying landslide susceptibility. The slope angle in the Rudraprayag district varies from 0 to 87° indicating gentle to very-steep gradient. A very steep slope intensifies the potentiality of a landslide (LS) event [40,41]. Duration and magnitude of sunray over the surface is not similar in all facets of slope, which in terms endures varying rate of physical weathering, vegetation growth and soil moisture thereby producing different potential segments for LS hazard. Therefore, slope aspect ( Figure 4c) is a vital LCF, particularly in the regions of higher elevation. In the case of curvature (Figure 4d), a rise of curvature value denounces the landslide probability. SPI (Figure 4e) in the process of landmass movement indicates the power of erosion by water which stimulates energy to disintegrate and transport the surface elements.
Hydrological LCFs selected for this study are rainfall, drainage density and topographical wetness index (TWI). For spatial mapping of the rainfall distribution, the precipitation data from 1901 to 2019 (119 years) was considered. The average rainfall map of the study area was prepared using the Kriging method based on the data collected from different stations. Rainfall (Figure 4f) is a widely used LCF in LSS mapping as it triggers the processes like soil slip, debris fall by declining the stability, integration and compaction of the surface forming materials. Towards the north and north-eastern part of the district, amount of average annual rainfall gradually increases from 150 cm to 300 cm. The arrangement of the drainage is the outcome of long-term interactions among geology, topography, climate and lithology. Generally, high intensity of drainages epitomizes low infiltration and prompt run-off and vice versa [52,53]. A large share of the inventory landslides was found along the high drainage density regions (Figure 4g). TWI (Figure 4h) is controlled by topographic factors and the index values represent the tendency of flow accumulation at a specific point and ability of the gravitational forces to flow the water downward in a slope segment [54].
Soil related LCFs in this research are soil texture and soil depth. Six categories of soil texture ( Figure 4j) were mapped i.e., coarse loamy, fine loamy mixed loamy, loamy skeletal, loamy sandy, and loamy. Loamy skeletal to the gravelly loamy texture of soil has more erosion susceptibility as it is found in the regions of higher relative relief and drainage density. Higher soil depth with higher elevation in the slope segment designates a higher possibility of the LS hazard (Figure 4i).
Geological factors of the present study are distance to lineaments, geology and earthquake zonation. Areas near to the fault lines are geologically high unstable compared to the distant locations. The lineaments were derived from analyzing satellite images using environment for visualizing images (ENVI) and PCI Geomatica software and from Bhuvan web platform of ISRO. Distances were categorized into five classes (Figure 4k) by applying Jenk's natural breaks classification scheme. The geological distinctions of this district are Chamoli Quartzite, Tourmaline Granite, Chandpur, Nagthat, Tourmaline Granite, Gneiss-Magmata, Granite 500 Ma, Chail-Ramgarh-Amri, and Shalkhalas-Jutogh (Figure 4l). Numerous past landslide events had occurred in Gneiss-Magmata and Chamoli Quartzite category regions. Occurrence of landslides is also influenced by the earthquake events, seismic disturbances, mountain slope direction and dip direction near the faults. Therefore, the earthquake zonation of the district was included as an LCF in this research. The zones were categorised according to MSK intensity scale (Medvedev-Sponheuer-Karnik scale) in Figure  4m. This scale (Table 2) helps in evaluating the gravity of the seismic event based on the effect on the district.
Human related factor i.e., road density was considered for the landslide susceptibility mapping. Construction of roadways on the steep slopes removes the support of above landmasses. Removal of the landmass across the slope segment interrupts the stability along the gradient resulting in topographic failure event. Therefore, road density is a vital factor in landslide susceptibility mapping. The landslide inventories showed a tendency of taking place alongside the roadways (Figure 4n).
Environmental LCFs considered in this study are Normalized Difference Vegetation Index (NDVI) and Land Use/Land Cover (LU/LC). NDVI (Figure 4o) is considered as it is the reflection of diverse character of soil properties i.e., moisture, structural composition, climatic event and vegetation types, which are allied with landmass removal processes [51]. Various LU/LC categories exhibit different role in LS hazard in the downhill such as agriculture area and settlement have a positive relation with the landslide (Figure 4p). The district was classified into several broad LU/LCs which are evergreen forest (47.31%), deciduous forest (8.53%), settlement (1.50%), cropland (5.68%), grazing area (6.53%), barren rocky land (8.76%), perennial water bodies (2.05%), glacier region (1.05%), seasonal water body (0.79%), and permanent snow region (15.27%).  and is an extended form of the model CART (Classification and Regression Tree). In regression and classification difficulties, RF model is recommended as a strong tool as illustrated in the literature [12,26,56]. The RF model is an ensemble one and modified bootstrap aggregation [57]. This model represents the modified bagging or bootstrap aggregating method and is therefore ensemble in nature. Integrated decision rules of RF model generates a very large number of trees to form a 'forest' and these are grown on the basis of bootstrapped sample applying CART framework with a random subgroup of variable chosen at every node. The final discernment about class membership and model building are decided on the basis of election of all decision trees [58]. In the current research, to run the RF model, 'RandomForest' package of 'R Studio' was used along with ArcGIS 10.3.1.

ANN
Based on the actual construction and behaviour of physiological neurons of the nervous system, ANN (artificial neuron networks), an artificially structured model computes process and transmits information to another level to provide response or output. The salient advantages of this model are: it differentiates and recognizes several sets of data within a large range of data, does not require existing experience, knowledge, or a pre-existing frame to train data [59,60]. In the present study, MLP (multi-layer perceptron) architecture of ANN was employed which have three layers, i.e., input, hidden, and output functioned through several neurons.
The MLP architecture in the present analysis operates in the following way: (a) assignment of random weight to all the linkages of causative landslide factors to start the model, (b) find out the activation rate using inputs and linkages, (c) find the error rate at the output node and recalibrate all the linkages between hidden nodes and output nodes, (d) cascade down the error to hidden nodes, (e) further recalibration and repetition of the process till convergence criterion is met, and (f) finally uses the final linkage weight score to the output node to produce the result. (2) where are the inputs, is corresponding weights and is the output derived through the function of .

SVM
Support vector machine is a binary classifier of supervised learning based on the law of structural risk minimization [61] in the field of data mining. This technique differentiates the hyperplane construction from training sub-set of inventory data. Under the original space of n coordinates (xi factor in vector x), differentiation of hyper-plane was produced between the points of two distinctive classes [24]. SVM constructs a classification hyper-plane in the middle of the highest margin as it revealed the highest limit of separation among the classes [62]. The classification is defined as the +1 (point over hyper-plane) which represent landslide presence pixels and -1 (point under the hyper-plane) designates landslide absence pixels. The training sub-sets which are adjacent to optimal hyper-plane are termed as support vectors. After the acquisition of a decision surface, the categorization of the new data can be prepared, including the training sub-set, i.e., label pairs ( ) with ∈ , ∈ 1, 1 and 1. . . . . . , . [63]. For producing the LSMs in the present study, X (vector space) was comprised of altitude, slope gradient, slope aspect, curvature, stream power index, rainfall, drainage density, topographical wetness index, distance from fault, geology, earthquake zonation, distance to roads, normalized difference vegetation index, and land use/land cover. The target of support vector machine is to extract optimal differentiating hyper-plane in which the training sub-set could split into absence and presence of landslides [−1, +1]. In this study, for the preparation of a landslide model using SVM, the radial basis function (RBF) kernel was used. More detailed regarding the SVM function are available in the literature of Roy et al. [64].

Logistic Regression (LR)
The logistic regression involves analysing a problem in which the outcome is estimated based on binary response of the dependent variable i.e., true or false and 0 or 1, specifically presence or absence of landslide event in this analysis which is predicted by one or a set of causative variables [65]. The relative contribution of a number (n) of independent variables (V1, V2,…..Vn) on a dependent (Y) is used to predict a logit transformation of probability based on the presence or absence of landslide. Though the LR model does not directly define the susceptibility, but the derived probability values could be helpful to draw an inference. The model was fitted using Equation (4).
where p represents the probability of landslide event Y (dependent variable), p/(1−p) is the likelihood ratio, C0 is intercept and C1, C2,…,Cn are the coefficients which were used to compute the contribution of the independent variables V1, V2,…,Vn. In this analysis, considering landslide event (training dataset) as dependent variable and sixteen causative parameters as independent variable, the spatial relationship was calculated using Equations (3) and (4). Similar to Akgun [66], in this study, using the probabilities, the possibility of occurrence of landslide was acquired using the coefficients from Equation (5).  where B is the logistic regression coefficient value, while EarthquakeB, LULCB, GeologyB, Soil textureB, and AspectB are the logistic regression coefficient values.

Ensemble of Models
Ensemble modeling is a method of combining the effects of different models into a unified embedded model to improve predictive capacity [67,68]. This methodology has raised the attention of researchers, in particular those familiar with the models of data mining and machine learning [63,68,69]. Using the weighted aggregation of an individual model, the ensemble models can be generated. The way these weights are measured, however, is complex. The present research uses a form of integration defined as a heterogeneous framework that uses mathematics of multiplication, division, addition, and subtraction where an improved equation was developed for the analysis and goes further [57]. Using Equation (5), the weighted mean expression was used to construct the ensemble models within the four parameters.
where EM is the resulted ensemble model, AUSRCi is the area under the success rate curve (AUSRC) value of the ith single model (Mi).

Validation Techniques
The present research includes two methods for evaluating the models, i.e., discrimination accuracy measures and reliability accuracy measures.

Discrimination Accuracy Measures
Area under the curve, precision, and accuracy were computed to trace the discrimination ability of the models under landslide and non-slide classes. Performance evaluation of the models by ROC (receiver operating characteristics) is a proficient and widely accepted method in the field of landslide susceptibility mapping [64]. The area under curve (AUC) validates the competence and accuracy of the employed models [70]. ROC method has been used in several disciplines including medical studies, physics and other branches, however, it has been most frequently applied for validating the spatial models of hazard predictions [64,68]. The ROC curve evaluates the outcome of classifier methods over the whole range of its cut-off points [71]. The curve represents the reliability in twodimensional form and there remain two types of error in the simulations which are chance of presence or chance of absence of any event [72]. In the X-axis, the curve plots 1-specificity, and Y-axis represents sensitivity. The sturdiness of prognostic models and their steadiness has also been assessed using precision, and efficiency which have computed following Equations (7) and (8). Values of precision and efficiency nearer to 1 also specify the better prediction ability of the models.
The precision, efficiency, and AUC in the present work computed following the equations below: Efficiency TP TN TP TN FP FN (8) where TP denotes true positive and TN is the true negative, which reflects a perfect classification of pixels under landslide presence and landslide absence, respectively. FP expressed as false positive refers to a quantity of landslide absence pixels that are incorrectly categorized under landslide present classes. FN is the false negative and signifies the number of pixels with landslide presence incorrectly categorized under landslide absence class. AUC value nearer to the 1 specifies a stronger performance of the models [73].

Reliability Accuracy Measures
The reliability of the models was assessed using two statistical techniques viz. mean absolute error (MAE) and root mean square error (RMSE). The MAE can be elaborate as the summation of difference between actual and predicted values divided by all total observations of a dataset. The MAE and RMSE have been illustrated in Equations (10) and (11).
where n is the sample size for training or testing sub-set, Lpredicted and Lactual are the predicted and actual values of landslides (dependent factor). RMSE and MAE values less than 0.5 designate a better prediction ability of the models and greater than 0.5 resembles less ability [74].

Model Prioritization Using Compound Factor
The compound factor (CF) method seeks to allot consecutive rank to the variables depending on their delegacy relevance aiming the goal. The average values of consigned ranks to the variables epitomize their relative priority [75,76].
1 (12) where 'R' is the variable rank and is the number of variables. To identify the best fit model for producing LSMs the CF-based prioritization was performed in terms of precision, efficiency, AUC, MAE, and RMSE values among the 15 models.

Result of Relief-F Analysis
Selection of the appropriate attributes or LCFs to enhance the model's capability is essential before the calibration and Relief-F in this regard accredited as an effective factor selection method [77]. This method was involved to compute the average merit (AM) of the LCFs and AM value greater than '0′ of LCFs denotes higher importance of the factor while negative or values <0 indicate very low or negative significance. The importance of the LCFs is shown in Table 3 in which they are ordered according to descending average merit.

LSMs by Individual Models
In the first phase of analysis, landslide susceptibility maps (LSMs) were produced applying the single-model based assessment method. Therefore, four different LSMs of the Rudraprayag district were produced by random forest, artificial neural network, support vector machine and logistic regression models individually. All the models were functioned in 'R studio' using R programming along with SPSS v17. For all the models, the susceptibility was categorised into five classes based on Jenk's natural breaks classification algorithm such as Very-high (VH), High (H), Moderate (M), Low (L) and Very-low (VL). For the classification of LSMs into VH, H, M, L, and VL, four techniques were tested, i.e., geometric interval, natural breaks, equal interval, and quantile. The calibration results provided by the Jenk's natural break are more prominent in distributing landslides as it has no bias and minimum deviation in intra-class and maximum deviation in inter-class which was also found in the work of Arabameri et al. [78].
RF model based LSM (Figure 5b) is categorised into 18.68%, 15.85%, 17.63%, 17.91%, and 29.93% area of the district under VH, H, M, L, and VL classes, respectively. Random forest was fabricated by the relative weight of the mean decrease in accuracy of the Gini index. The ANN model produced 21.59%, 15.08%, 18.93%, 18.29%, and 26.11% areas are susceptible under VH, H, M, L, VL categories, accordingly (Figure 5a). Similarly, the LR model-based result categorized 23.62% area of the district as very highly susceptible, 13.31% as highly susceptible and 15.15%, 18.52%, and 29.38% areas as moderate, low, and very-low category, respectively (Figure 5c). In the case of the SVM model-based LSM (Figure 5d), the very-low (25.93%) and moderate (16.81%) susceptibility class cover the largest areas of the district respectively followed by very-high (23.83%), high (14.38%) and low susceptibility (25.93%) classes. The comparative visualization of the individual model based LSMs are represented in Figure 6.

LSMs by Two Ensemble Models
The second phase of the analysis was to produce LSMs using an ensemble of two-models (

LSMs by Ensemble of Three-Models
In the third phase, the ensemble approach was executed considering three models differently

LSMs by Ensemble of Four-Models
Finally, all four models, i.e., RF, ANN, LR, and SVM models, were combined to prepare the LSM of the Rudraprayag district. The ensemble model of ANN-LR-RF-SVM classified 21.52% area under the very high susceptibility category, 14.63% area under high, 16.47% area under moderate, 19.72% area under low, and 27.65% area of the district under the very low susceptibility category (Figure 9).

Results of the Validation Techniques
The validation method in this research includes two approaches such as reliability and discrimination ability evaluation. AUCROC, precision and efficiency were computed for measuring the discrimination abilities of the models while the reliability was assessed applying MAE and RMSE which were applied in various research works as validation techniques [79][80][81]. The present research using effective model-based approaches has produced fifteen LSMs considering both single and combined techniques. Therefore, precision, efficiency, AUC, MAE, and RMSE of 15 models using both of the training and testing sub-set GLIs were computed and presented in Tables 4 and 5

Result of Variable Importance Analysis
The results of the assessment of the variable importance based on the mean decrease of Gini coefficient using RF model are shown in Figure 11. In addition, the coefficients of logistic regression which also indicates the importance of the variables are shown in Table 6. As Figure 11 shows, all the causative factors generally contributed to the landslide susceptibility model. However, altitude is the most important factor, followed by drainage density, road density, lineament density and slope respectively. Like the mean decrease of Gini, the coefficients of logistic regression (Table 6) also indicate a strong relationship between the landslide occurrence and the altitude, rainfall, drainage density, road density and distance from lineament. These results showed that most of the landslides occurred in elevated areas with high drainage density and road density.

Discussion
Determining the LCFs, generation of LSMs and selection of the best-fitted model is the initial stage of landslide hazard assessment and the present research was intended to attempt that. As the landslide is the uncertainty of spatial association among the several factors, different models can be applied for the prediction. During the last decades, various empirical and numerical models have been attempted to forecast geohazards [14,79,82]. However, those models have some limitations and assumptions [15,24]. Recently, data mining models have been largely popularised in the research community, especially in modelling different environmental risks due to their ability to analyse complex connection between predictors (LCFs) and responses (landslide or non-landslide). Besides this, the ensemble of different data mining models along with statistical models also has been used effectively to produce LSMs [13,63,83]. This work was carried out through four phases of ensemble approaches that made the work novel compared to other works. In resolving the over-fitting problem, an ensemble of these machine learning models has more capability because these models automate the task of selecting multiple datasets without pre-assumption and can handle massive data [77]. Furthermore, the evaluation of the prediction ability of the models is a fundamental phase in any modelling work and it could be judged using both threshold dependent and independent methods [14,83]. Therefore, the AUC, precision, efficiency, MAE, and RMSE of each modelling approach were assessed considering both the training sub-set GLIs and testing or validation sub-set GLIs. Measurement of AUC including training sub-set represents success rate while the prediction rate was computed using testing sub-set of GLIs [14,73].
Considering the training sub-set GLIs, in the individual modelling, the AUC was found highest (85.12%) in ANN model followed by RF (84.59%), SVM (83.39%) and LR (82.17%) model (Table 4). In terms of precision and efficiency values, ANN leads the highest performance followed by SVM, RF and LR model (Table 4). However, both the reliability measures MAE and RMSE showed best results for RF model followed by ANN, SVM and LR model (Table 4). LR-RF showed the best result in terms of precision and efficiency when two models were combined and all the two-model ensembles have shown greater accuracy than the individual model in term of AUC ( Figure 10). When the three models were combined, ANN-LR-RF revealed the best result under each measure of precision (0.871), efficiency (0.878), AUC (87.83%), MAE (0.016) and RMSE (0.019) ( Table 4). At the last phase, to produce LSM, all four models were combined and the performances of the ensemble ANN-RF-SVM-LR model were: precision = 0.784, efficiency = 0.804, AUC = 87.76%, MAE = 0.151, and RMSE = 0.195. Correspondingly, the performances (reliability and discrimination ability) of the models were assessed considering testing sub-set of GLIs as well. Individually, the LR model regarded as the best in terms of precision (0.825) and efficiency (0.829) but AUC defined the RF model as the best fit for predicting landslide susceptibility. The RF model has greater acceptability in spatial risk mapping as found in various research works [30,54,78]. Similarly, the validation results of two-model ensemble, three-model ensemble and four-model ensemble have been listed in Table 5.
Compound factor-based prioritization of models in terms of the discrimination and reliability accuracy measures considering training and testing datasets was used to find the best-fitted models. Regarding the success rate of the produced LSMs, the CF method showed that an ensemble of ANN-LR-RF models was the best fit followed by RF-SVM-LR and ANN-RF-SVM-LR models, respectively. In the context of prediction performance, priority was the highest for ANN-LR-RF again followed by ANN-RF-SVM-LR and RF-SVM-LR models. The produced LSMs by the individual and ensemble models were classified into five susceptibility classes. The areas covering different susceptibility classes were varied from model to model. The best fit ensemble model, i.e., ANN-LR-RF showed that 15.74%, 12.67%, 19.33%, 21.97%, and 30.28% areas have VH, H, M, L, and VL susceptibility for the landslide. However, the output of the work will help the planners for the management of landslide in this area.

Conclusions
The key field of concern in this study is the analysis of landslide susceptibility. Single and ensemble models were successfully used for landslide susceptibility evaluation in the Rudraprayag district of Garhwal Himalaya. ROC curve, efficiency, accuracy, MAE, and RMSE were employed to assess and compare these LSS models. Finally, using the CF the best model was selected. The findings of the validation techniques show that the LSS models prepared by ANN, FR, SVM, LR, and their ensembles performed well in this analysis. The ensemble of ANN, RF, SVM, and LR models played an improving role as they improved the values of AUC. However, with respect to the results of robustness, the ANN-RF-LR ensemble had the highest reliability in terms of results where the highest agreement was found in the efficiency, accuracy, AUC of success rate curve, and predictive rate curve values, MAE and RMSE. The present research shows that ensemble models are effective tools that could be used in landslide prediction and landslide management for efficient decision-making by the authorities. In the present analysis, these models were utilized in the Rudraprayag district. However, for wider use, the output of such models may be used for areas of similar geo-environmental conditions. The key drawback of this analysis is that we used a limited set of LCFs and a fixed training testing samples dataset ratio (70:30) for the all the models. The performance of the models will have to be further tested through a variety of training and testing sets of sample datasets as well as LCF sets.