Policy Research Working Paper 10264 Machine Learning and Sensitivity Analysis Approach to Quantify Uncertainty in Landslide Susceptibility Mapping Mohammed Basheer Thomas Oommen Masatsugu Takamatsu Sachi Suzuki Urban, Disaster Risk Management, Resilience and Land Global Practice December 2022 Policy Research Working Paper 10264 Abstract Mitigating the impacts of landslides requires quantifying Tracts and Sylhet divisions of Bangladesh to understand the the susceptibility of different infrastructures to this hazard implications of weight uncertainty for road susceptibility through landslide susceptibility mapping. The mapping to landslides. The case study results show that distance to requires overlaying the spatial effects of multiple factors roads is the most influential factor to determine the like- that contribute to the occurrence of landslide events (rain- lihood of the occurrence of landslide events, followed by fall, land cover, distance to roads, lithology, and slope) the land cover type. Given weight uncertainty, the percent- and this process requires assigning weights to the different age of road lengths in the study area under extremely high factors contributing to landslides. This study introduces a susceptibility to landslides ranges from around 20 to 38 per- new statistical approach for quantifying the weights used cent. The tolerance level to weight uncertainty is a crucial in landslide susceptibility mapping and their associated determinant of investment costs and is ultimately a critical uncertainty. The proposed approach combines machine element for decision making to relevant institutions and learning (random forest classification) with large-scale sen- affected stakeholders. A conservative selection of weights sitivity analysis to derive the uncertainty ranges of weights from within the uncertainty range (a weight combination used in landslide susceptibility mapping. The study demon- that results in the highest susceptibility) means that the risk strates the approach for a case study of the Chittagong Hill is minimized but with a high investment cost. This paper is a product of the Urban, Disaster Risk Management, Resilience and Land Global Practice. It is part of a larger effort by the World Bank to provide open access to its research and make a contribution to development policy discussions around the world. Policy Research Working Papers are also posted on the Web at http://www.worldbank.org/prwp. The authors may be contacted at mohammedadamabbaker@gmail.com, toommen@mtu.edu, mtakamatsu@worldbank.org, ssuzuki1@worldbank.org The Policy Research Working Paper Series disseminates the findings of work in progress to encourage the exchange of ideas about development issues. An objective of the series is to get the findings out quickly, even if the presentations are less than fully polished. The papers carry the names of the authors and should be cited accordingly. The findings, interpretations, and conclusions expressed in this paper are entirely those of the authors. They do not necessarily represent the views of the International Bank for Reconstruction and Development/World Bank and its affiliated organizations, or those of the Executive Directors of the World Bank or the governments they represent. Produced by the Research Support Team Machine Learning and Sensitivity Analysis Approach to Quantify Uncertainty in Landslide Susceptibility Mapping Mohammed Basheer, Thomas Oommen, Masatsugu Takamatsu, Sachi Suzuki JEL Classification Q54 Climate • Natural Disasters and Their Management • Global Warming Keywords Landslide Risk Reduction, Road Asset Management, Disaster Risk Management, Geographic Information Systems, Disaster Risk Assessment 1 1. Introduction Landslides are a natural process of Earth’s geological formation and occur when gravitational force overcomes the resistance force of land materials (Turner, 2018). However, when landslides occur near human settlements and economic activity areas, they can result in losses of lives and livelihoods. Landslides are expected to increase due to the impact of human activities, such as land use or land cover modifications, hill cutting, and climate change (Crozier, 2010). Landslides can be initiated by several factors, such as rainfall, seismic activity, volcanic eruption, snowmelt, stream erosion, coastal erosion, changes in surface water levels, changes in groundwater, and/or human activities (Schuster and Wieczorek, 2002). The Global Landslide Catalog (GLC) reports around 6,800 rainfall-triggered landslide events since 2007, with around 20% of these recorded landslide events located in the South Asia Region (Kirschbaum et al., 2010). GLC reports over 25,000 fatalities that resulted from rainfall- triggered landslides since 2007. Landslides result in direct and indirect costs to humans (Turner, 2018). The direct costs include the loss of lives, properties, and infrastructure, while the indirect costs are associated interruptions to economic activities (Turner, 2018). Landslide susceptibility mapping is an approach used to provide a spatial classification of potential landslide areas (Fell et al., 2008; Sarkar and Kanungo, 2004). Landslide susceptibility maps are often produced using a heuristic method called “qualitative map combination,” in which the effects of multiple factors that contribute to landslide occurrence are combined using weights (Fell et al., 2008). These weights are based on expert knowledge and are challenging to determine due to varying expert opinions. Such a method has an unpredictable error margin, as individual perceptions can influence expert choices. Behavioral scientists have highlighted notable biases in decisions involving financial or life losses (Tversky 2 and Kahneman, 1981). Thus, the fuzzy logic approach has been applied in landslide susceptibility mapping to address the uncertainty in expert opinion on factors contributing to landslide occurrence (Feizizadeh et al., 2014; Tangestani, 2009). Although the fuzzy logic helps characterize and incorporate the overall uncertainty in expert opinion, it does not eliminate or reduce possible expert biases. To address this issue, statistical approaches have been used to establish statistical relationships between the factors that contribute to landslides and the occurrence of the landslide phenomena (Ayalew and Yamagishi, 2005; Gómez and Kavzoglu, 2005; Lee et al., 2004; Youssef et al., 2016). Lee et al. (2004) used an artificial neural network to determine the weights of the factors contributing to landslide occurrence. While statistical methods have contributed to reducing the uncertainty in landslide susceptibility mapping stemming from subjective choices, they have not accounted for nor quantified the uncertainty in the statistical determination of weights due to landslide event sampling for training and testing the statistical models. In 2019, the World Bank commissioned a study on geohazard risk management and resilient road asset management in Bangladesh. The study, referred to as “the Phase 1 study” for the rest of this paper, was conducted by IMC Worldwide Ltd —an international development consultancy— and the Centre for Environmental and Geographic Information Services (CEGIS). The study aimed to assess the landslide risk of roads in the Sylhet and Chittagong Hill Tracts (CHT) divisions of Bangladesh, and to prepare mitigation and adaptation investment plans for the 20 most critical landslide hotspots. Based on expert judgment, the study used a qualitative approach to overlay the different factors contributing to rainfall- triggered and earthquake-triggered landslides in the study regions. Table 1 reports the weights used by the Phase 1 study for conducting rainfall-triggered landslides susceptibility analysis. The Phase 1 study included reviewing landslide risk literature, conducting landslide risk 3 assessment, conducting a vulnerability assessment, developing resiliency plans, and providing training to stakeholders. The relevant stakeholders in Bangladesh involved in the Phase 1 study were the Local Government Engineering Department (LGED), the Roads and Highway Department (RHD), and the City Corporation of Chattogram and Sylhet. Table 1 The weights used by IMC Worldwide and CEGIS for landslide susceptibility mapping of rainfall-triggered landslides. Factor Weight Distance to Road 0.14 Rainfall 0.14 Land Use 0.40 Lithology 0.24 Slope 0.08 Building on the data collected during the above-described Phase 1 study, this paper proposes a statistical approach for determining the weights used for combining the input factors for landslide susceptibility mapping, alongside the uncertainty resulting from sampling landslide events used in weight quantification. We used the same data, factor classes, and scores of the Phase 1 study to be able to directly evaluate the outcomes of using our proposed statistical approach for determining weights compared to the qualitative approach used in the previous study. Our proposed approach combines a machine learning model (Random Forest Classification) and large-scale sensitivity analysis for sampling the testing and training data of the machine learning model. Data on factors contributing to landslide susceptibility are used to train and test machine learning models for predicting past landslides. Sensitivity analysis is conducted to sample the locations of past landslides used for training the models. The relative importance of factors contributing to landslide susceptibility is determined for each sensitivity iteration, with the relative importance of all sensitivity experiments representing the weight uncertainty range. We demonstrate the proposed approach for a case study of road susceptibility 4 to rainfall-triggered landslides in the CHT and Sylhet divisions of Bangladesh. Bangladesh is located in South Asia and faces a high risk of landslides. Notably, the CHT and Sylhet divisions are more vulnerable to landslides than other regions in Bangladesh, with many recent rainfall- triggered landslides (Ahmed, 2015). Results show that considering the uncertainty of weight estimation in landslide susceptibility analysis enables probabilistic assessment of the landslide hazard and can help select an uncertainty tolerance level for decision-making. Case study results show that, on average, distance to existing roads is the highest contributor to landslide occurrence, followed by the land cover type. Given the uncertainty in weight estimation, the percentage length of local roads and highways under high or extremely high susceptibility can range from around 15% to 35%. 2. Study area The study area is Bangladesh's CHT and Sylhet divisions (Fig. 1). Bangladesh is located in South Asia and has a land area of around 130,000 km2 and a population of around 165 million as of 2020 (The World Bank, 2022). India borders the country from the west, north, and east, and Myanmar in the southeast, while the Bay of Bengal borders the south. Bangladesh is administratively divided into eight divisions. CHT and Sylhet divisions are located in the eastern part of Bangladesh and occupy around 25% of the country's land surface area. The CHT Division is characterized by rugged hilly topography, with many hogback ridges. The weather in the CHT Division is tropical, with mean annual rainfall ranging from around 2,500 mm to approximately 3,800 mm. The Sylhet Division has a less rugged topography compared to the CHT Division (see Fig. 1) and has an average annual rainfall of around 3,300 mm. In the last two decades, many landslides have occurred in Bangladesh, causing loss of lives and livelihoods and economic, environmental, and infrastructural damage (Ahmed, 2021). 5 Rabby and Li (2019) developed an inventory of landslides for the Chittagong Hilly Areas of Bangladesh and reported 730 landslide events from 2001 to 2017, with 356 events in the summer of 2017. The number of landslide events in Bangladesh has been increasing in recent years due to increasing rainfall intensity as a result of climate change (Khan et al., 2020), deforestation (Reddy et al., 2016), urbanization, and hill cutting (Ahmed, 2021). Fig. 1 Location and topography of the Chittagong Hill Tracts (CHT) and Sylhet divisions in Bangladesh. The transportation infrastructure in Bangladesh includes networks of local roads and highways. The RHD of Bangladesh is responsible for constructing and maintaining major highways and bridges, whereas LGED is responsible for local-level infrastructure. The CHT and Sylhet divisions include networks of roads and highways of around 19,300 km and 14,600 km, respectively. This road network is vulnerable to natural hazards such as landslides. 6 Interruption to road services due to landslides can cause direct economic losses due to infrastructure damage and indirect losses that result from the impact on socio-economic activities that rely on roads and highways. To address the landslide issue in Chittagong, Islam et al. (2017) and Ahmed et al. (2018) developed landslide susceptibility mappings. Islam et al. (2017) demonstrated that quality mapping can be made through open-source spatial data, multi-criteria evaluation, and experts’ knowledge. Ahmed et al. (2018) compared four methodologies to develop landslide susceptibility maps, including (a) machine learning methods – the random forest, artificial neural network, (b) support vector machine, (c) multiple linear regression, and (d) principal component analysis and multiple linear regression, and concluded (a) as the best model for developing the landslide early warnings, validating by receiver operating characteristics. This study will use machine learning methods to quantify the weights that are used for landslide susceptibility mapping that informs decision making of investments for protection. 3. Methodology Fig. 2 shows a flowchart of the hybrid machine learning and sensitivity analysis approach introduced in this study for determining weight uncertainty in landslide susceptibility mapping. The approach consists of seven interlinked, iterative steps. These steps are numbered from 1 to 7 in Fig. 2. 3.1. Trigger and factor data Spatially distributed data for the factors and triggers contributing to landslide occurrence are collected in the first step of the proposed approach. The factors considered in this study 7 (which focuses on rainfall-triggered landslides) are rainfall, land cover, distance to roads, lithology, slope, and distance to road. Fig. 3 shows maps of the five factors. Spatially distributed data of the distance to the nearest road for each grid within the study area was created (Fig. 3a) based on vector data of roads and highways in the CHT and Sylhet divisions. The road and highway vector data were acquired from RHD and LGED authorities of Bangladesh. A land cover map for the year 2015 (Fig. 3b) was obtained from CEGIS. CEGIS also provided gridded geological data for the study area (Fig. 3c). A slope layer (Fig. 3d) was created based on a 30-meter Digital Elevation Model (DEM) obtained from the NASA Shuttle Radar Topography Mission (SRTM; Jarvis et al., 2008). Spatially distributed annual rainfall data with a 30-year return period (Fig. 3e) were generated based on gage observations from Bangladesh Meteorological Department. The gaging stations considered for generating the rainfall layer are Chittagong, Comilla, Cox'sbazar, Feni, Kutubdia, Rangamati, Sandwip, Shrimangal, Sitakunda, Sylhet, and Teknaf. 8 Fig. 2 Hybrid machine learning and sensitivity analysis approach for landslide susceptibility mapping. The layers of the five factors were resampled and aligned at a 50-meter spatial resolution. The layers were then reclassified based on scores estimated by CEGIS for a previous study. Table 2 reports the scoring scheme used to reclassify the layers of the five factors that contribute to landslide occurrence. These classes and their scores where developed by IMC Worldwide Ltd and the Centre for Environmental and Geographic Information Services (CEGIS) during a study commissioned by the World Bank. For land cover, the locally-obtained land use data from 9 2015 was used for broad classifications into the 5 categories: Agriculture, Water bodies, Forest, Bare soil based on judgement by the local consultant firm, CEGIS. For example, only natural (or reserved) forest was categorized into the Forest category, while agriculture, orchards, vegetation, and other plantation were categorized into the Agriculture category as vegetation areas are mostly used to homestead or cultured vegetation for food and livelihood ran by local communities. The shifting cultivation was categorized under the Human Interventions in addition to roads and settlements. We acknowledge that these classes and scores could be refined, but we used them to be able to compare our results with the results of this previous study. 10 Fig. 3 Factors contributing to landslide susceptibility considered for the CHT and Sylhet divisions. Table 2 Score assignment for different factor classes for the CHT and Sylhet divisions. Score Distance to road (km) Land cover Lithology Slope (degree) Rainfall (mm) 1 >1,000 Bare soil and St. Martin Limestone 0-5 0-600 water bodies 11 Score Distance to road (km) Land cover Lithology Slope (degree) Rainfall (mm) 2 750-1,000 Forest Valley Alluvium and 5-10 600-700 Colluvium, Marsh Clay and Peat, Estuarine deposit 3 500-750 N. A. N. A. 10-15 700-750 4 N. A. N. A. N. A. N. A. 750-800 5 N. A. N. A. Girujan Clay Formation, N. A. N. A. Tidal Mud Deposit, Bokabil 6 250-500 Agriculture N. A. 15-20 800-850 7 N. A. N. A. N. A. N. A. 850-900 8 N. A. N. A. N. A. N. A. N. A. 9 0-250 Human Dihing, Dupitila and Tipam >20 >900 intervention 3.2. Landslide and non-landslide locations inventory The second step of the proposed approach is to compile an inventory of past landslide locations and develop an inventory of non-landslide locations (i.e., locations where no landslides occurred in the past). Point-based data for past landslides in the CHT and Sylhet divisions were compiled from multiple sources, including Ahmed et al. (2014), Comprehensive Disaster Management Programme II (2012), and GLC. The green points in Fig. 4 show the locations of past landslides in the inventory. Random and spatially diverse non-landslide points were created in the CHT and Sylhet divisions such that: a) The number of non-landslide points is approximately equal to the number of landslide points in each division. A roughly equal number of landslide and non-landslide points were ensured to avoid bias in training the machine learning model, as explained later. b) Each non-landslide point is at least 200 m away from the nearest landslide location. 12 Fig. 4 Inventory of landslide and non-landslide locations for the CHT and Sylhet divisions. Google Earth satellite images were used to ensure the randomly generated non-landslide points are not located at past landslide locations that are not captured in the landslide inventory. If a point was suspected to be located over a past landslide location based on satellite images, the point was removed from the non-landslide locations inventory. The red points in Fig. 4 show the non-landslide locations included in the inventory. The inventory includes 301 13 landslide locations, of which 282 are in the CHT Division and 19 are in the Sylhet Division. A total of 283 non-landslide locations are included in the inventory. A total of 268 non-landslide locations are in the CHT Division, while 15 are in the Sylhet Division. 3.3. Sampling of locations and extraction of factors In the third step of the proposed approach, the score values (see Table 2) of the factors considered are extracted at the landslide and non-landslide locations compiled in the second step. The extracted data were prepared in a tabular format with multiple columns. One of the columns includes binary values (zero or one) representing whether a location is a landslide or a non-landslide. The rest of the columns show the score value for each factor. The total number of rows equals the total number of landslide and non-landslide locations compiled in the inventory. In the fourth step of the proposed approach (Fig. 2), a fraction of the landslide and non- landslide locations (compiled in the second step) is randomly selected. The rationale behind this step is to reflect the fact that not all past landslide and non-landslide locations are included in the inventory. The use of a limited inventory of landslide and non-landslide locations for landslide susceptibility analysis is a common issue in Bangladesh and many regions worldwide. Many landslides occur in remote areas and are often unrecorded on the ground. Basing statistical quantification of weights for landslide susceptibility mapping on a limited fraction of landslide and non-landslide locations brings uncertainty that needs to be considered in risk analysis and infrastructure investment. The proposed approach enables repeating this random selection step of locations multiple times through sensitivity analysis to quantify the uncertainty in weights stemming from limited observed data. In this analysis, we randomly selected a sample of 80% of landslide and non-landslide locations in each sensitivity iteration, as shown 14 in Fig. 2. We conducted sensitivity analysis around the splitting percentage locations for training. We test 60%, 70%, and 90% in addition to the selected 80% value. Table 3 shows the implications for splitting assumptions for weight uncertainty. Table 3 1Sensitivity of factor weights to splitting of testing and training data. Factor Road distance Rainfall Land cover Lithology Slope Data split 90/10 80/20 70/30 60/40 90/10 80/20 70/30 60/40 90/10 80/20 70/30 60/40 90/10 80/20 70/30 60/40 90/10 80/20 70/30 60/40 Max 0.29 0.32 0.33 0.34 0.24 0.27 0.27 0.28 0.28 0.29 0.30 0.34 0.20 0.21 0.22 0.23 0.16 0.18 0.20 0.21 weight Min 0.21 0.20 0.17 0.17 0.18 0.16 0.15 0.14 0.19 0.17 0.15 0.14 0.14 0.13 0.12 0.12 0.12 0.11 0.10 0.1 weight Median 0.25 0.25 0.25 0.25 0.21 0.21 0.21 0.20 0.23 0.22 0.22 0.22 0.17 0.17 0.17 0.18 0.14 0.14 0.15 0.15 weight 3.4. Machine learning and factor weights In the fifth step of the proposed approach, we used machine learning to calculate the weights of the factors that have been identified in the first step in determining the occurrence of landslides. We used the Random Forest Classification Machine Learning Algorithm (Breiman, 2001) to develop classification models for predicting whether a location is a landslide or a non-landslide location. To this end, we used the data sampled in the fourth step to train a Random Forest Classification model in each iteration. The model target (i.e., dependent variable) is to predict a binary value (zero or one) representing whether a location is a landslide or a non-landslide location. The model features (i.e., independent variables) are the score values of the factors identified in the first step. A random forest model is an ensemble of decision trees; each tree is created based on a random selection from the training data and the available features (Svetnik et al., 2003). During model training, each tree in the forest is designed by selecting leaf nodes and splits to maximize prediction quality. The prediction of a random forest model is made by aggregating the predictions of individual trees in the forest based on the majority vote (i.e., the mean) (Breiman, 2001). In this study, we used the Gini criterion to measure the 15 quality of splits and guide model training. We also used 100 tree predictors (i.e., individual trees in the forest), each with a maximum depth of 100 nodes. Once a Random Forest Classification model is trained based on the data sampled in the fourth step of the proposed approach, model performance is tested on the remaining data that have not been used in training the model. The Overall Accuracy metric (Oommen et al., 2010) was used to evaluate the performance of the model with the training and testing data. To calculate the Overall Accuracy, confusion matrices were created based on model performance with the training and testing data. A confusion matrix is a table that counts the number of True Positives (TP), True Negatives (TN), False Positives (FP), and False Negatives (FN) in model predictions. Equation 1 shows how the overall accuracy can be calculated based on a confusion matrix. The Overall Accuracy value can range from 0 to 1, with 0 indicating that none of the model predictions is true, while a value of 1 indicates that all model predictions are true. Thus, high Overall Accuracy values are desirable. + Overall Accuracy = +++ (1) The Overall Accuracy performance of the model produced in the fifth step is evaluated against a performance threshold. This threshold sets the overall acceptable level of accuracy in model predictions. For this case study application, we used an overall accuracy threshold of 0.80. As Fig. 2 shows, if model performance with either the testing or training data (as measured by the Overall Accuracy) is below the specified threshold, the model is discarded, and the process is repeated from the fourth step, where a new random sample of landslide and non- landslide locations is drawn. If the model performance is satisfactory, we proceed to the sixth step of the proposed approach. 16 In the sixth step, the importance of features for model predictions is calculated for the Random Forest Classification model developed in the fifth step. Feature importance reflects the weight of factors in determining the occurrence of a landslide event. Feature importance is a metric that can be calculated for Random Forest Classification models. For each model feature (i.e., independent variables), the feature importance metric can take any value from 0 to 1, but the summation of the importance values for all features equals 1. The higher the value of the importance metric of a feature, the stronger that feature is in influencing predictions. A feature importance value of 1 indicates that the feature is the only one contributing to model predictions, while a feature importance value of 0 means that the feature does not influence model predictions. The importance values of the landslide factors represent the weights by which each factor contributes to the occurrence of a landslide event. In the Bangladesh case study, feature importance was calculated based on impurity. This approach performs well in our case as the features considered have a low cardinality. After the feature importance values are calculated, they are stored before proceeding to the next iteration. Multiple iterations of steps 4-6 of the proposed approach are performed to capture the uncertainty in determining factor weights resulting from using a limited inventory of landslide and non-landslide locations, as explained earlier. A target number of iterations should be selected. After the sixth step, if the target number of iterations is reached, the analysis proceeds to the seventh step; if not, a new iteration starts by returning to the fourth step (see Fig. 2). In our case study application, we specified a target number of iterations of 10,000. After the target number of iterations is achieved, the range within which the weight of each factor varies (based on different iterations) is used to represent the uncertainty range. The uncertainty range of the weights of the different factors is used in the seventh step to derive uncertainty in landslide susceptibility mapping, as explained in the next section. 17 3.5. Landslide susceptibility mapping under factor weight uncertainty In the last (seventh) step of the proposed approach, sensitivity analysis is conducted around road susceptibility to landslides. Weight values are randomly selected within the weight range computed in the sixth step of the framework, with all values within the range equally probable. We conducted 80 random weight draws for this case study to cover the weight uncertainty range determined in the sixth step. The randomly selected weights are used to overlay the reclassified factor layers generated in the first step. The resulting composite layer represents the Landslide Susceptibility Index (LSI). Equation 2 describes how the reclassified factor layers are overlaid. = ∑ =1 × (2) Where LSI is the Landslide Susceptibility Index, is the weight of the factor number i, is the score layer of the factor number i, and n is the number of factors considered. The resulting LSI layers were classified into very low, low, moderate, high, and extremely high classes. Table 4 reports the LSI classes used in this study. These classes were obtained from CEGIS; CEGIS created this classification scheme for a previous landslide susceptibility assessment in Bangladesh. After classifying the susceptibility layer, the length of local roads and highways passing through each landslide susceptibility class is calculated in addition to the percentage lengths of roads and highways within each class. Table 4 Landslide Susceptibility Index classes adopted in this study. Landslide susceptibility Landslide Susceptibility Index class (LSI) Very low 0 ≤ LSI ≤ 2 Low 2 < LSI ≤ 4 Moderate 4 < LSI ≤ 5 High 5 < LSI ≤ 6 18 Landslide susceptibility Landslide Susceptibility Index class (LSI) Extremely high 6 < LSI 4. Results and discussion 4.1. Machine learning and weight factor uncertainty For the CHT and Sylhet divisions, Fig. 5a shows the performance of the random forest classification models with the training and testing data. The figure shows that the Overall Accuracy of some machine learning models reached as high as 0.91 with the training data and 0.90 with the testing data. Model performance varies due to random sampling of landslide and non-landslide locations used for training and testing, as explained in 3.4 and Fig. 2. Fig. 5 Random forest classification models for weight estimation for the CHT and Sylhet divisions: (a) performance of the models; (b) box plots of the uncertainty in the weights of factors contributing to landslide susceptibility. The ends of the boxes in panel b represent the upper and lower quartiles, the continuous vertical lines inside the boxes mark the medians, the circles show the outlier points, and the whiskers extend to the maximum and minimum values, excluding the outliers. Overall, Fig. 5a shows the machine learning models perform better with the training data than with the testing data. The lowest recorded Overall Accuracy with the training data was around 0.87, while the lowest Overall Accuracy with the testing data was lower at 19 approximately 0.80. Generally, the performance with the training and testing data is approximately similar and ranges between 0.8 and 0.91, indicating the model does not overfit the training data. Fig. 5b shows boxplots of the uncertainty range of the weights of the factors contributing to rainfall-triggered landslides in the CHT and Sylhet divisions of Bangladesh. These uncertainty ranges resulted from the iterations over steps 4-6 of the proposed approach (see Fig. 2), as explained in 3.4. The results in Fig. 5b show that distance to roads is the highest contributor to landslide events compared to the other factors considered. The weight of distance to roads took a range from around 0.20 to roughly 0.30, with a median of approximately 0.25. This strong influence of distance to roads can be explained by the fact that around 75% of the past landslides included in the inventory occurred within 250 m of a road. Without proper planning and design, hill cutting for road construction can reduce slope stability and contribute to the occurrence of landslide events (Chisty, 2014). The hybrid machine learning and sensitivity analysis approach revealed that land cover type is the second strongest influencer on the occurrence of landslide events, with a weight that ranges between 0.18 and 0.28 with a median weight of around 0.18. Around 50% of the landslide events included in the inventory occurred in areas with a land cover class of “human interventions” followed by “agricultural land,” in which around 41% of landslide events occurred. As Fig. 5b shows, slope is the weakest factor influencing the occurrence of landslide events in the CHT and Sylhet divisions. The weight by which slope contributes to landslide events ranges from around 0.12 to approximately 0.18, with a median of around 0.14. This weak influence can be explained by the low number of past landslide events recorded in areas with high slopes. In the landslide inventory, only 2 landslide events occurred in areas with slopes greater than 20 degrees, whereas 228 events (75% of events in the inventory) occurred in areas with less than 10 degrees slopes. 20 4.2. Road susceptibility to landslides under factor weight uncertainty Fig. 6 shows a parallel coordinate plot of the sensitivity analysis of factor weight values within the uncertainty range determined using machine learning. A parallel coordinates plot is a visualization approach to present multi-dimensional data (Inselberg, 2009). Fig. 6 shows how varying the weight values within their uncertainty ranges impacts the proportion of road lengths in the CHT and Sylhet divisions that fall under high or extremely high susceptibility to landslides. Fig. 6 Parallel coordinates plot of weight combinations within the uncertainty range and their resulting road susceptibility to landslides in the CHT and Sylhet divisions. Each grey line in Fig. 6 connects the weight values (axes in the black box on the right side) and their corresponding road susceptibility status (axes in the black box on the left side). The thick red and blue lines mark weight combinations with the highest and lowest road susceptibility to landslides, respectively. Results show that given weight uncertainty, the 21 percentage of road lengths under extremely high susceptibility to landslides can range from around 20% to 38%. Higher uncertainty was revealed in the percentage of roads under high susceptibility to landslides, with an estimate ranging from around 28% to 38%. The primarily crossing grey lines between the two axes in the black box on the left side of Fig. 6 indicate a weak correlation between the percentage of road lengths under high and extremely high susceptibility to landslides, given weight uncertainty. Also, the random changes in line directions between all axes of Fig. 6 indicate the non-linear nature of the interplay between the factors contributing to the occurrence of landslides and road susceptibility to landslides. This non-linearity highlights the important role that machine learning and sensitivity analysis can play in understanding weight uncertainties in landslide susceptibility mapping. The thick red line in Fig. 6 shows that the highest estimated extremely high road susceptibility to landslides results from weight values of 0.18, 0.24, 0.28, 0.17, and 0.13 for rainfall, land cover, distance to roads, lithology, and slope, respectively. Contrarily, the thick blue line shows that a weight combination of 0.24, 0.18, 0.23, 0.18, and 0.17 for rainfall, land cover, distance to roads, lithology, and slope, respectively, results in the lowest estimate of the percentage of road length in the CHT and Sylhet divisions under extremely high susceptibility. The weight values shown by the thick red and blue lines in Fig. 6 confirm the non-linear nature of the interplay between weights and the resulting road susceptibility estimate. For instance, although distance to road has been revealed by machine learning to be the most influential factor on landslide occurrence, reducing the weight of this factor to the lowest value within the uncertainty range does not necessarily result in the lowest susceptibility. The lowest susceptibility estimate also involves increasing the weights of less influential factors like rainfall and slope. The spatial distribution of the factor values and roads and highways shape this non-linear complexity. 22 4.3. Decision-making under uncertainty regarding road susceptibility to landslides Numerical simulation models can have uncertainties related to their parameterization and/or structures (Reed et al., 2022). One of the goals of this study is to present an approach that can estimate a reasonable range of parameter weightage used for risk assessment simulations, which would help practitioners and decision makers to make informed investment decisions acknowledging the uncertainties. As noted by Ullman (2001), robust decision-making accounts for the uncertainties. However, the level of uncertainty that can be tolerated would impact investment and management decisions (Regan et al., 2005). Fig. 7 shows the percent exceedance of the fraction of road lengths in the CHT and Sylhet divisions under extremely high susceptibility to landslides. This figure was produced based on the 80-weight sensitivity iteration conducted in the seventh step of the proposed approach (Figure 2). The 80 sensitivity experiments were performed by randomly drawing weight combinations from within the uncertainty range, assuming all values are equally probable. Fig. 7 shows, dependent on the set of weight combination selected, total length of “extremely high risk” road varies from 20% to 37%, and the resultant risk maps look different (Figure 8). In one simulation, only around 20% of the length of roads and highways would fall under extremely high susceptibility to landslides, potentially resulting in low investment cost to mitigate the risk of landslides. On the other hand, the length of roads and highways under extremely high susceptibility to landslides increases to around 37%, suggesting a higher investment cost would be required. All other simulations fall between the two extreme scenarios. Practitioners or decision makers can select a suitable weight combination depending on the overall risk involved (i.e., after factoring in vulnerability) and the resources available to implement mitigation measures. This selection process should involve different related governmental institutions and representatives of stakeholders affected by road susceptibility to landslides (e.g., civil society 23 and the private sector). A conservative investment plan would allow no tolerance for uncertainty (i.e., red point in Fig. 7 or red line in Fig. 6). This would ensure a robust road infrastructure with respect to landslides but would require a high investment budget. Fig. 7 Percent exceedance of the fraction of road length under extremely high susceptibility in the CHT and Sylhet divisions given weight uncertainty. The point colors (i.e., red and blue points) correspond to the lines with similar colors in Fig. 6. Fig. 7 can be used in two different ways. First, if a hypothetical $100 million investment budget has been allocated, Fig. 7 can be used to assess the level of tolerance and robustness that can result from this investment budget. For instance, if the allocated budget is enough to invest in 25% of roads under extremely high susceptibility (on the x-axis of Fig. 7), this translates into around 60% uncertainty tolerance (on the y-axis of Fig. 7). The second use case for Fig. 7 is to determine the required investment budget based on a preselected level of tolerance. To further unpack the implications of selecting a tolerance level for weight uncertainty, Figure 8 shows the spatial distribution of road susceptibility to landslides in the CHT and Sylhet divisions with the two extreme examples explained above (i.e., full and no tolerance levels). The figure shows that the spatial distribution of roads and highways under high or extremely 24 high susceptibility to landslides changes depending on the adopted tolerance level. The higher the level of tolerance, the lower the length of roads and highways under high or extremely high susceptibility. Although this study demonstrated the methodology to identify road susceptibility, accuracy of the result is limited due to the limitation of available data. Lack of comprehensive landslide data resulted in the reduced relative importance of slope towards landslide occurrence. It is reflected as the lowest weight given to slope in the figure 6. To improve this model and support other similar studies, improvement in data collection is necessary. Data on landslides that occur in remote steep areas are needed to improve the accuracy in estimating the weights of the factors that contribute to the occurrence of landslide events. In addition to data on the locations of past landslides, other data related to the impacts of past landslides (e.g., extent, number of damaged houses, and number of affected people) can be used to train machine learning models to predict vulnerability to landslides. 5. Conclusions Weight uncertainty in mapping susceptibility to landslides can have substantial investment implications. Quantifying and understanding this uncertainty is a key in selecting an acceptable uncertainty tolerance level for robust decision-making. In this study, we propose a new approach that combines machine learning with large-scale sensitivity analysis in an attempt to better quantify factor weights in landslide susceptibility mapping and their associated uncertainties. The proposed approach is demonstrated for a case study of road susceptibility to rainfall-triggered landslides in the CHT and Sylhet divisions of Bangladesh. The case study results reveal that distance from existing roads and highways is the most influential factor in the occurrence of landslide events, highlighting the vital role that 25 geotechnical and engineering design of hill cutting can play in reducing the landslide hazard. Given the uncertainty in factor weights, the percentage length of roads and highways in CHT and Sylhet divisions that fall under extremely high susceptibility can range from approximately 20% to 38%. Selecting an acceptable level of tolerance would enable adopting a road susceptibility level within this uncertainty range, which can be used to inform risk analysis and investment decisions. 26 Fig. 8 Susceptibility of roads and highways in the CHT and Sylhet divisions based on combinations of weights that result in (a-b) the lowest susceptibility and (c-d) the highest susceptibility. The frame colors of the maps (i.e., red and blue) correspond to the lines with similar colors in Fig. 6. This study focused on one source of uncertainty: the weights used to overlay the factors contributing to landslide events. However, other potential uncertainty sources should be investigated. For instance, in the case study application, a scoring scheme for different factor classes has been adopted from a previous study conducted by CEGIS (see Table 2). However, CEGIS scores are based on expert judgment, which can result in biases. Another potential source of uncertainty in landslide susceptibility mapping is the susceptibility classes, which were obtained from CEGIS and can introduce uncertainty. Sensitivity analysis of landslide susceptibility mapping to factor scores and hazard classes should be conducted alongside factor weights to fully account for the uncertainties involved. In the case study application, we used factor data resampled at a 50-meter spatial resolution. Using a finer spatial resolution could improve performance and reduce the uncertainty around the weights used in landslide susceptibility mapping. For instance, improving the spatial resolution of the digital elevation model layer would improve slope estimates and could influence the resulting weights of the slope layer. In this study, we relied on an inventory of landslide locations compiled by IMC Worldwide and CEGIS. While this inventory comprehensively accounts for many past landslide locations, other locations in remote areas are typically not recorded due to inaccessibility or the absence of infrastructure at risk. Any bias in the landslide inventory would propagate to a bias in the weight uncertainty analysis. For example, if most past landslide locations were recorded in urban areas, the weight of the land cover layer would increase, as areas with human interventions are assigned a high score in reclassifying the land cover layer. Improving the landslide inventory can improve weight estimates and reduce biases. 27 The weight uncertainty results produced for Bangladesh show marked differences from the weights adopted in the previous study conducted by IMC Worldwide and CEGIS. These differences are because our analysis used a statistical approach while the previous study used a qualitative approach based on expert judgment. This difference highlights a potential bias from solely relying on expert judgment in weight estimation. Ideally, a weight estimation procedure combining statistical and qualitative approaches, whereby weight uncertainty is determined using the statistical approach (i.e., machine learning and large-scale sensitivity analysis), then the qualitative approach (i.e., expert judgment) is used to select a weight combination (or multiple combinations) from within the uncertainty range. Such a combination of statistical and qualitative approaches could further reduce the uncertainty range and provide more guidance to decision-making on road investment. The next step is to translate the uncertainty quantification and transform it to the decision- making support tool. Although accounting for different uncertainty sources can help make robust investment decisions, understanding and interpreting large-scale sensitivity analysis results can be challenging for decision-makers who need to choose an acceptable level of uncertainty tolerance. New data visualization techniques such as parallel coordinates plots (Inselberg, 2009) can help decision-makers and stakeholders sift through the multi-dimensional uncertainties related to landslide investments. Online data visualization tools such as Polyvis (https://www.polyvis.org/) can support inter-institutional and multi-stakeholder deliberation on the uncertainties related to landslide susceptibility and help zero in on an acceptable level of uncertainty tolerance.* * An interactive parallel coordinates plot of the sensitivity analysis results shown in Fig 6 can be accessed on Polyvis via the following link, where different weight combinations can be filtered, and their susceptibility implications can be explored: https://www.polyvis.org/sheet/pvphofbsoj8f. 28 6. Code availability statement The Random Forest Classification machine learning model is open-source and was used through the Python package scikit-learn: https://github.com/scikit-learn/scikit-learn. Sensitivity analysis was performed using the Python application programming interface of the Gdal package: https://gdal.org/. 7. References Ahmed, B., 2021. The root causes of landslide vulnerability in Bangladesh. Landslides 18, 1707–1720. https://doi.org/10.1007/s10346-020-01606-0 Ahmed, B., Rahman, M.S., Islam, R., Sammonds, P., Zhou, C., Uddin, K. and Al-Hussaini, T.M., 2018. Developing a dynamic Web-GIS based landslide early warning system for the Chittagong Metropolitan Area, Bangladesh. ISPRS International Journal of Geo- Information, 7(12), p.485. Ahmed, B., 2015. Landslide susceptibility mapping using multi-criteria evaluation techniques in Chittagong Metropolitan Area, Bangladesh. Landslides 12, 1077–1095. https://doi.org/10.1007/s10346-014-0521-x Ahmed, B., Rahman, M.S., Rahman, S., Hug, F.F., Ara, S., 2014. Landslide Inventory Report of Chittagong Metropolitan Area, Bangladesh. Ayalew, L., Yamagishi, H., 2005. The application of GIS-based logistic regression for landslide susceptibility mapping in the Kakuda-Yahiko Mountains, Central Japan. Geomorphology 65, 15–31. https://doi.org/10.1016/j.geomorph.2004.06.010 Breiman, L., 2001. Random Forests. Machine Learning 45, 5–32. https://doi.org/10.1023/A:1010933404324 Chisty, K.U., 2014. Landslide in Chittagong City: A perspective on hill cutting. Journal of Bangladesh Institute of Planners 7, 1–17. Comprehensive Disaster Management Programme II, 2012. Landslide Inventory & Land-use Mapping, DEM Preparation, Precipitation Threshold Value & Establishment of Early Warning Devices. Crozier, M.J., 2010. Deciphering the effect of climate change on landslide activity: A review. Geomorphology 124, 260–267. https://doi.org/10.1016/j.geomorph.2010.04.009 Feizizadeh, B., Shadman Roodposhti, M., Jankowski, P., Blaschke, T., 2014. A GIS-based extended fuzzy multi-criteria evaluation for landslide susceptibility mapping. Computers and Geosciences 73, 208–221. https://doi.org/10.1016/j.cageo.2014.08.001 Fell, R., Corominas, J., Bonnard, C., Cascini, L., Leroi, E., Savage, W.Z., 2008. Guidelines for 29 landslide susceptibility, hazard and risk zoning for land-use planning. Engineering Geology 102, 99–111. https://doi.org/10.1016/j.enggeo.2008.03.014 Gómez, H., Kavzoglu, T., 2005. Assessment of shallow landslide susceptibility using artificial neural networks in Jabonosa River Basin, Venezuela. Engineering Geology 78, 11–27. https://doi.org/10.1016/j.enggeo.2004.10.004 Inselberg, A., 2009. Parallel coordinates: Interactive visualisation for high dimensions, in: Liere, R., Adriaansen, T., Zudilova-Seinstra, E. (Eds.), Trends in Interactive Visualization: State-of-the-Art Survey. Springer London, London, pp. 49–78. https://doi.org/10.1007/978-1-84800-269-2_3 Islam, M.A., Murshed, S., Kabir, S.M., Farazi, A.H., Gazi, M.Y., Jahan, I. and Akhter, S.H., 2017. Utilization of open source spatial data for landslide susceptibility mapping at Chittagong District of Bangladesh—an appraisal for disaster risk reduction and mitigation approach. International Journal of Geosciences, 8(04), p.577. Jarvis, A., Reuter, H., Nelson, A., Guevara, E., 2008. Hole-filled SRTM for the globe Version 4 [WWW Document]. CGIAR-CSI SRTM 90m Database. URL http://srtm.csi.cgiar.org (accessed 12.30.18). Khan, M.J.U., Islam, A.K.M.S., Bala, S.K., Islam, G.M.T., 2020. Changes in climate extremes over Bangladesh at 1.5 °C, 2 °C, and 4 °C of global warming with high-resolution regional climate modeling. Theoretical and Applied Climatology 140, 1451–1466. https://doi.org/10.1007/s00704-020-03164-w Kirschbaum, D.B., Adler, R., Hong, Y., Hill, S., Lerner-Lam, A., 2010. A global landslide catalog for hazard applications: Method, results, and limitations. Natural Hazards 52, 561– 575. https://doi.org/10.1007/s11069-009-9401-4 Lee, S., Ryu, J.H., Won, J.S., Park, H.J., 2004. Determination and application of the weights for landslide susceptibility mapping using an artificial neural network. Engineering Geology 71, 289–302. https://doi.org/10.1016/S0013-7952(03)00142-X Oommen, T., Baise, L.G., Vogel, R., 2010. Validation and Application of Empirical Liquefaction Models. Journal of Geotechnical and Geoenvironmental Engineering 136, 1618–1633. https://doi.org/10.1061/(asce)gt.1943-5606.0000395 Rabby, Y.W., Li, Y., 2019. An integrated approach to map landslides in Chittagong Hilly Areas, Bangladesh, using Google Earth and field mapping. Landslides 16, 633–645. https://doi.org/10.1007/s10346-018-1107-9 Reddy, C.S., Pasha, S.V., Jha, C.S., Diwakar, P.G., Dadhwal, V.K., 2016. Development of national database on long-term deforestation (1930-2014) in Bangladesh. Global and Planetary Change 139, 173–182. https://doi.org/10.1016/j.gloplacha.2016.02.003 Reed, P.M., Hadjimichael, A., Malek, K., Karimi, T., Vernon, C.R., Srikrishnan, V., Gupta, R.S., Gold, D.F., Lee, B., Keller, K., Thurber, T.B., Rice, J.S., 2022. Addressing uncertainty in multisector dynamics research. https://doi.org/10.5281/zenodo.6515285 Regan, H.M., Ben-Haim, Y., Langford, B., Wilson, W.G., Lundberg, P., Andelman, S.J., Burgman, M.A., 2005. Robust decision-making under severe uncertainty for conservation management. Ecological Applications 15, 1471–1477. https://doi.org/10.1890/03-5419 30 Sarkar, S., Kanungo, D.P., 2004. An Integrated Approach for Landslide Susceptibility Mapping Using Remote Sensing and GIS 70, 617–625. Schuster, R.L., Wieczorek, G.F., 2002. Landslide triggers and types, in: Proceedings of the 1st European Conference on Landslides, Prague. Balkema, Rotterdam. pp. 59–78. Svetnik, V., Liaw, A., Tong, C., Christopher Culberson, J., Sheridan, R.P., Feuston, B.P., 2003. Random Forest: A classification and regression tool for compound classification and QSAR modeling. Journal of Chemical Information and Computer Sciences 43, 1947– 1958. https://doi.org/10.1021/ci034160g Tangestani, M.H., 2009. A comparative study of Dempster-Shafer and fuzzy models for landslide susceptibility mapping using a GIS: An experience from Zagros Mountains, SW Iran. Journal of Asian Earth Sciences 35, 66–73. https://doi.org/10.1016/j.jseaes.2009.01.002 The World Bank, 2022. The World Bank Database [WWW Document]. URL https://data.worldbank.org/ (accessed 5.20.21). Turner, A.K., 2018. Social and environmental impacts of landslides. Innovative Infrastructure Solutions 3, 1–25. https://doi.org/10.1007/s41062-018-0175-y Tversky, A., Kahneman, D., 1981. The framing of decisions and the psychology of choice. Science 211, 453–458. Ullman, D.G., 2001. Robust decision-making for engineering design. Journal of Engineering Design 12, 3–13. https://doi.org/10.1080/09544820010031580 Youssef, A.M., Pourghasemi, H.R., Pourtaghi, Z.S., Al-Katheeri, M.M., 2016. Landslide susceptibility mapping using random forest, boosted regression tree, classification and regression tree, and general linear models and comparison of their performance at Wadi Tayyah Basin, Asir Region, Saudi Arabia. Landslides 13, 839–856. https://doi.org/10.1007/s10346-015-0614-1 8. Acknowledgments The authors would like to thank GFDRR-EU which funded the South Asia Region Building Resilience to Landslide and Geohazard Risk in the South Asia Region Program and Africa Fellowship Program for funding the work. 31