Multi-criteria risk evaluation by integrating an analytical network process approach into GIS-based sensitivity and uncertainty analyses

ABSTRACT Geographic information system (GIS)-based multi-criteria decision analysis (MCDA) is commonly used to solve a range of complex spatial problems. We have incorporated an analytical network process (ANP) as a central part of MCDA and used this approach for land subsidence susceptibility mapping (LSSM). We used an ANP approach in combination with the concepts of global sensitivity and uncertainty analyses, in order to minimize any associated errors. This approach was tested for a LSSM application in north-western Iran, using a methodology that consisted of three distinct phases. An ANP approach was used in the first phase to derive criteria weightings by analysing relevant criteria, which included topographic, hydrologic, climatologic, geological, and anthropological indicators. The second phase involved evaluating the uncertainty and sensitivity of areas susceptible to subsidence as a function of the derived weightings, using Monte Carlo simulations. The third phase then used a land subsidence inventory database to validate the results and to measure any improvement in accuracy following sensitivity analysis. Results indicated a 6% improvement in the accuracy of the ANP method as a result of using an integrated global sensitivity analysis.


Introduction
Geographic information system (GIS)-based multi-criteria decision analysis (MCDA) is a powerful modelling methodology in the spatial sciences. GIS and MCDA have been successfully combined to solve a variety of spatial problems (Feizizadeh and Ghorbanzadeh 2017). The GIS-MCDA method combines data derived from a number of geographical factors into a single index of evaluation (Chen et al. 2010a). GIS-MCDA also forms a framework in which to manage the different components of any spatial decision problem, and can be used to evaluate relationships between the different elements in complex problems (Malczewski 2006;Feizizadeh et al. 2014a).
The analytic network process (ANP) is one of the most common GIS-MCDA techniques and has been widely used to derive criteria weightings in modelling tasks (Saaty 1996a). This method is a variation of the analytic hierarchy process (AHP) that was developed in the context of network approaches in order to minimize the errors involved in a traditional AHP (Chen and Yang 2011). The ANP is considered advantageous over other MCDA methods because decision-makers, especially those dealing with spatial problems, are able to observe the interdependencies in a 'real world' problem and assess reliability through clusters and a network structure (Taslicali and Ercan 2006). A number of previous publications have described the use of ANPs in spatial decision-making problems (Neaupane and Piantanakulchai 2006;Chen and Yang 2011;Ferretti 2011). However, this approach requires decision-makers with highly relevant expertise in order for them to be able to generate pairwise comparison matrices in a professional and effective manner (Taslicali and Ercan 2006). The results of an ANP are highly dependent on decisions and assessments provided by the experts, which can be a major source of uncertainty in such analytical network decision-making processes.
Decision makers in GIS-MCDA procedures make use of a variety of spatial data, take into account the opinions of experts, analyse the data, and set preferences according to rules previously established for the decision-making process (Malczewski 2004;Rahman et al. 2012). GIS-MCDA has been widely used in a variety of modelling tasks such as, for example, landslide susceptibility mapping, land suitability assessment, site selection, land subsidence susceptibility mapping, etc. In this research, we have used GIS-MCDA to determine the weightings for criteria relevant to land subsidence susceptibility, for use in LSSM.
Although they yield accurate results in most cases, it has been demonstrated that there is a degree of uncertainty associated with GIS-MCDA methods, which can sometimes lead to inaccurate results (Feizizadeh et al. 2014a). GIS-MCDA results have therefore not always been completely satisfactory for decision makers, because of this inherent uncertainty Feizizadeh and Kienberger 2017). Criteria weightings can contribute significantly to the uncertainty in GIS-MCDA methods because the user's individual preferences with regard to the input data are assumed to be actual rules in the spatial decision-making process, even though some of these preferences may in fact be erroneous. Furthermore, errors in the input data and the level of uncertainty in the structure of the modelling are not entirely independent of each other and can interact in a number of different ways (Loucks et al. 2005). The exact sources and magnitudes of uncertainties in the results obtained from a GIS-MCDA model are therefore typically unclear and difficult to compute, due to the heterogeneity of input data and the variety of parameters involved . Ignoring this uncertainty can have an unfavourable influence on the accuracy of the results and any maps that are subsequently produced (Sarmento et al. 2008). A number of attempts have been made to improve the accuracy of MCDA as it is potentially a valuable tool for solving a wide range of spatial problems (e.g. Boroushaki and Malczewski 2008;Rahman et al. 2012).
The authenticity of spatial simulation models and GIS-MCDA techniques depends on how the relative importance of each parameter has been calculated. Without a sensitivity analysis and a good understanding of the model's structural integrity and the accuracy of the input factors, the accuracy of any prediction made by the model will typically remain uncertain (Norton 2015).
Despite recent progress in GIS-MCDA, there is still no single function available for assessing the level of uncertainty in the results. There is therefore a great need for a method to analyse and evaluate these errors and to assess the reliability of the results obtained . In response to this need, concepts of sensitivity and uncertainty have been developed that aim to determine the sources of error in GIS-MCDA in order that they can then be minimized. An uncertainty analysis is not the same as a sensitivity analysis and a number of researchers have integrated the two into GIS-MCDA (Crosetto et al. 2000;Ligmann-Zielinska and Jankowski 2012;Feizizadeh et al. 2014b;Norton 2015;Erlacher et al. 2017). The two terms (uncertainty and sensitivity) refer to different concepts: an uncertainty analysis attempts to describe the entire set of possible outcomes and their associated probabilities of occurrence, while a sensitivity analysis attempts to determine the magnitude of change in model output that results from a minor change in model input. A sensitivity analysis therefore measures the change in model output with respect to a localized region of specific inputs. The same set of model runs can, however, often be used for both uncertainty and sensitivity analyses (Loucks et al. 2005).
Sensitivity analyses investigate the magnitude of a model's response to different sources of variation, such as in the input data-set, the parameters selected, and the assumptions made (Tenerelli and Carver 2012). Such analyses allow a level of confidence to be determined in the model output, which helps decision makers to recognize and measure confidence intervals in the modelling results . One of the main advantages of uncertainty and sensitivity analyses in GIS-MCDA is their ability to estimate confidence intervals by evaluating the responses visible in the modelling results to small changes in the input data (Crosetto et al. 2000;Feizizadeh et al. 2014b). They are also able to reveal the relationships between input and output parameters in a multi-criteria decision problem (Chen et al. 2010b). The importance of uncertainty and sensitivity analyses in validating the outcomes of GIS-MCDA is being increasingly acknowledged (Chen et al. 2010a) and should be considered an essential part of any GIS-MCDA application.
The authors have previously discussed the importance of sensitivity and uncertainty analyses in GIS-MCDA and demonstrated that the variety of possible errors associated with criteria weighting based on AHP can be minimized through such analyses, thereby significantly improving the accuracy of results Feizizadeh et al. 2014b). Building on the results of this earlier work we have aimed to address the uncertainties associated with particular components of GIS-MCDA (such as criteria weighting through ANP), and to define a unified approach to uncertainty and sensitivity analyses by means of a global sensitivity analysis (GSA). We have therefore applied sensitivity and uncertainty analyses to an ANP method and determined the various reactions in order to minimize uncertainty and improve accuracy. The ANP method was used to take advantage of its criteria-weighting capabilities, as will be discussed in Section 3.2. In the next section, we review the uncertainties associated with ANP results. The subsequent sections then describe the implementation of spatially explicit sensitivity and uncertainty analysis (SESUA) methods and the generation of LSSM products on the basis of both classic approaches and integrated approaches. The results are then validated through a case study on land subsidence susceptibility in Marand County, Iran.

Study area and data-set
The case study area was Marand County in the East Azerbaijan province of north-western Iran ( Figure 1). Marand County has an area of 3,286 km 2 and is important for its natural landscape and as a population centre. The altitude reaches more than 3250 m above mean sea level in mountainous areas to the north and south, decreasing towards the central areas down to a low of just under 910 m in the river beds. The county is characterized by a range of landscape and soil types that make it suitable for a variety of agricultural activities, including both irrigated and rain-fed agriculture and orchard cultivation, which provide the main sources of income for the rural population (Khorrami 2016).
The study area is potentially vulnerable to land subsidence. The Regional Water Organization (RWO) reports that this region is one of six critical areas in the East Azerbaijan province in which any exploitation of groundwater (GW) resources is prohibited due to a dramatic fall in GW levels over recent years (Fakhri et al. 2015;Ghorbanzadeh 2016). Land subsidence is a common problem on most of the plains of Iran, and especially the north-western part of the country where agriculture is largely dependent on GW rather than on rainfall. The widespread use of GW resources within the study area makes it prone to land subsidence: it can take quite a few years from the instigation of land subsidence before damage becomes visible at the surface (Lee et al. 2012). Land subsidence can occur for a number of reasons and there can sometimes be a combination of triggers, leading to severe consequences. This has already occurred within our study area, where the massive exploitation of GW over the last three decades was a key factor, along with others, in the triggering of a land subsidence event (Fakhri et al. 2015).
For our research, we used a data-set comprising topographic, hydrologic, geological, anthropological, and climatic factors. We grouped these factors into clusters comprising a number of criteria related to active land subsidence that were selected on the basis of previous investigations. Eleven causal clusters were ultimately selected for use in this study, based on the following five main factors: & Climatic factors such as precipitation & Geological factors such as distance to fault and lithology & Hydrologic factors such as streams and depth of GW & Anthropological factors such as distance to aqueduct, land use, and distance to wells & Topographic factors such as land capability, slope, and DEM (see Table 1)  We grouped the factors into related clusters because clusters play an important role in the ANP method, as is explained in Section 3.2. Following all the necessary geometric and thematic editing, topologic corrections and simple basic calculations (e.g. buffers around faults and streams) were made and then all data-sets arranged in ArcGIS software as a raster layer with a resolution of 30 m, for further spatial analysis.

Methods: workflow
We first investigated the use of uncertainty analysis to calculate the variability of results, followed by a sensitivity analysis that subdivided the variability between the various types of input data. An uncertainty analysis first helps to reduce the uncertainties in how an MCDA operates and then parameterizes the stability of its outputs. This is achieved quite simply by marginally changing the input data and then evaluating the results produced by these changes, as previously discussed by Feizizadeh et al. (2014b). These investigations aimed to evaluate the sensitivity and uncertainty of GIS-MCDA for LSSM by (a) using sensitivity analysis to assess the accuracy of the MCDA technique, and (b) validating the certainty of outputs considering the improved accuracy of MCDA by applying sensitivity and uncertainty analysis. The methodology involved three main steps: criteria weighting/ ranking, sensitivity and uncertainty analyses, and validation ( Figure 2), as described below.
In the first step, all of the specified criteria (see Section 2 and Table 1) were assessed using the ANP method in order to derive appropriate weightings for the aggregation required to produce an LSSM map. This 'classic approach' to incorporating an ANP into GIS-MCDA is based on a multi-criteria evaluation that identifies those areas susceptible to subsidence on the basis of eleven possible causal factors. The second step focused on an 'integrated approach' by applying sensitivity and uncertainty analyses to the weightings derived using an ANP. Error propagation was simulated using the SESUA method, and a Monte Carlo simulation (MCS) was used to obtain the uncertainty associated with the results of the ANP (see Section 3.4). MCS is a frequently used practical method for evaluating uncertainty and variability (Firestone et al. 1997) and is suitable for use when no precise result can be computed by modelling (Singh and Arora 2015). This step was followed by the final step (described below), in which a variance-based GSA was used to determine the uncertainties in the weightings.
In the final step, LSSM maps were produced based on the weightings derived from both classic and integrated approaches. The results obtained using each approach were validated against a land subsidence inventory database using the relative operating characteristics (ROC) method. Figure 2 depicts the research methodology.

Training data-set and standardization
In order to determine the spatial correlations between each of the criteria and the decision-making of the SESUA model, it is necessary to use a number of data points with a normal distribution within the study area. For this objective, 23 GPS locations for known land subsidence were randomly selected from the 46 GPS locations available in the land subsidence inventory data-set. The remaining 23 locations were used for the validation. This data-set was actually based on initial signs of land subsidence, which were recorded in an extensive field survey over the study area. These data points served as input training data in the SESUA approach, for which they were attributed with eleven criteria that were prepared as separate GIS layers. Each of these criteria was represented on a map to be used for further analysis in the LSSM decision model. After attributing the training data points with each of the criteria, the values obtained were standardized to render all criteria on the same scale as a fuzzy standard. In other words all map cells were assigned a value between 0 and 1 on the basis of either benefit or cost factors. Those map cells that are more susceptible to land subsidence were assigned high standardized values (close to 1) and less susceptible cells were assigned low values (Rahman et al. 2012).

ANP: structure and function
The ANP method was introduced by Saaty (1996a) as an effective way of reducing the weaknesses and errors in a traditional AHP used for decision-making (Chen and Yang 2011). Even though the AHP method has been criticized by many researchers for certain deficiencies (Feizizadeh and Ghorbanzadeh 2017), it has still been used for criteria weighting and ranking in many GIS-MCDA applications. The ANP was developed as a complementary technique to the traditional AHP, with the former delivering accurate results more often than the latter (Gashti and Shahrodi).
Since the structure of an AHP is that of a hierarchical, top-down decision model, the criteria are assumed to be independent of each other and problems can arise if some of the criteria are in fact influenced by each other. The ANP model was developed to deal with this condition of interrelated elements in a network system (Pourebrahim et al. 2010, Mahdavi et al. 2013). According to Saaty (1996b), an ANP provides a logical way of dealing with dependence; it can therefore solve more complex problems than other MCDA multi-criteria decision-making models with interdependent relationships. As a general rule, the ANP technique can be subdivided into three different stages which are as follows (Zabihi et al. 2015): (1) Selecting the most suitable alternative in situations with conflicting and interrelated criteria by modelling the decision problem.
(2) Prioritizing the interrelated criteria through a framework, in order to evaluate it.
(3) Assessing the indirect influences of different criteria and their relative importance to one another through the network.
Since the ANP is considered to be a mathematical theory, it can methodically overcome all types of dependencies. The first step generates a number of pairwise comparison matrices of the criteria. The relative significance value can be determined on a scale of 1-9 (shown in Table 2), ranging from equal significance of variables to extreme dominance of one variable over another (Zabihi et al. 2015;Saaty, 1996b). The ANP forms a network that includes small to large matrices of criteria and possible alternatives (both called elements), grouped into clusters. The elements can be interlinked in any number of possible ways within the network. This network includes feedback and interdependence connections, within and between clusters (Aragon es-Beltr an et al. 2008). (1) The priority of each pairwise comparison matrix is calculated separately. The results of all individual pairwise comparison matrices are collected into an overall super matrix (1), in order to assess the final weightings of the elements (Saaty 1996b). The super matrix W illustrates the influence of each element on other elements within the network system (Aragon es-Beltr an 2008). By applying the ANP technique, we aim to calculate criteria weightings (overall priorities) from a limited matrix and assigning them on the basis of the effectiveness of the criteria. The limited matrix was obtained by calculating the square root of each weighted supermatrix by multiplying it by itself, in order to capture the overall priorities of the criteria. The outputs of the limited matrix therefore demonstrate the relative importance and hence the weighting of each criteria. The general form of the limited super matrix has previously been broadly described by Saaty (1996b). Within this complex process, the ANP enables more accurate modelling than other methods of complex problems, such as spatial problems. For any reciprocal pairwise comparison matrix, this approach requires (a) the consistency to be tested, and (b) the detection and modification of inconsistent element(s) (Ergu et al. 2014).
Although there are a number of different methods for evaluating consistency, the consistency ratio (CR) is the most popular and practical consistency index to use for pairwise comparison matrices in an ANP (Saaty 1990). The CR is expressed as where λ max is the maximum eigenvalue of the matrix A; RI (see Table 3) is the random index, which is related to the number of criteria; (n) is the order of the matrix A that is being compared with the pairwise comparison matrix (Malczewski and Rinner 2015). If CR < 0:1 the matrix has an acceptable level of consistency. As mentioned in Section 2, the case study for our research involved LSSM. In order to calculate the criteria weightings using an ANP, questionnaires were designed that could be used to obtain the expert knowledge required for the selection of ranking criteria. A number of experts in the fields of geology, geography, and hydrology were asked to complete the questionnaires and a total of 40 responses received that could then be used to determine criteria ranking through expert knowledge. The weightings ascribed to the 11 criteria were determined using the ANP technique in Super Decision software. Judgment matrices were structured for all clusters and the weightings of each matrix computed separately. The cluster matrix and the final weighted, unweighted, and limited matrices were established based on the relative importance of the criteria (see Tables 4-7).
In order to calculate the unweighted matrix (Table 5), pairwise comparisons were made between all criteria in a particular cluster that influenced each other and were dependent on each other. The unweighted matrix therefore involved all of the weighting vectors that were calculated by pairwise comparisons.   Anthropological Land use 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 Slope 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 0.14 The weighted matrix was produced by multiplying, element by element, the weightings from the unweighted matrix by the weightings of the related cluster (Neaupane and Piantanakulchai 2006). In the following, the weighted matrix was reached to limiting power so that the elements of the matrix to be convergent.
Calculation of the limited matrix has been previously explained in detail by Saaty (2008). The limited matrix (Table 7) shows that the 'distance to well' criterion had the highest weighting (0.21). As mentioned previously (in the description of the study area), agriculture is the most important activity in Marand County. However, agriculture in this area is heavily dependent on the wells from which GW is extracted. The 'distance to well' criterion is considered to be the most important criterion with regard to land subsidence because of the extreme lowering of the water table caused by over-extraction of GW. The subsequent reduction in pressure within the aquifer leads to an increase in the effective stress between different underground layers (Chai et al. 2004). The next most important criteria after the 'distance to well' were the 'land use' and 'slope' criteria. The CR value for the pairwise matrix in our study was calculated to be 0.053.

Uncertainty and sensitivity in ANP weightings
One approach to uncertainty analysis that is applicable to MCDA involved obtaining a sense of the error or uncertainty in the predictions from the uncertainty in the criterion weightings (Benke and Pelizaro 2010). The uncertainty associated with criteria weightings in a multi-criteria problem is thus based on the judgment of experts or stakeholders . Uncertainty is very common in spatial problems because of their complexity, as well as for a number of other reasons such as the ranking of factors on the basis of Eigen values derived from MCDA, but the absence of probabilities for each alternative can result in confusion for decision makers (Hsu and Pan 2009;Feizizadeh et al. 2014b). Since the judgments of a number of different experts are used in decisionmaking processes and their opinions will, in most cases, vary considerably, uncertainty is always likely to be present in the process and hence also in the subsequent results of MCDA, even if all judgments are considered to be correct. As discussed further in Section 4 below, the ANP method is more reliable than other MCDA methods (such as AHP) for criteria weighting/ranking in complex decision-making problems, i.e. in identifying the importance of each criterion. However, since the ANP relies on expert opinions, this method is assumed to be subject to a certain amount of bias in the pairwise comparisons of the criteria, and any inaccuracies can be easily conveyed into the weightings (Kritikos and Davies 2011;Feizizadeh and Blaschke 2012). Benke et al. (2008) listed some of the previous studies in which the uncertainties of various models have been addressed using different regressions and probabilistic methods, such as Markov Chain Monte Carlo (MCMC) methods (Kuczera and Parent 1998), or MCMC/Bayesian inference approaches (e.g. MCMC global sensitivity analysis (MCMC-GSA); Kanso 2006) and a Bayesian approach to total error analysis (BATEA; Kavetski et al. 2006). In order to deal with the sources of uncertainty, previous researchers (e.g. Hsu and Pan 2009;Benke and Pelizaro 2010;Feizizadeh and Blaschke 2012) have proposed integrating MCS into GIS-MCDA methods. The integration of MCS with MCDA involves using a probabilistic characterization of pairwise comparisons to tackle uncertainty (Bemmaor and Wagner 2000;Hahn 2003). This approach has been described by Feizizadeh et al. (2014b) and the reader is referred to this publication for further information.

Implementation of ANP-MCS
Our approach was to use a MCS to analyse the uncertainty associated with weightings derived using an ANP technique. Since the uncertainty in the input data and the final outputs (weightings) can be characterized by probability distributions, it was possible to use the ANP-MCS approach to analyse the uncertainties. There are a number of statistical sampling approaches available for analysing uncertainties associated with MCDA models, of which MCS is one of the most commonly used . Such simulation is one of the most suitable ways of studying the propagation of uncertainty of spatial models, in the absence of a large amount of information concerning the functional form of the inaccuracies (Eastman 2003;Tenerelli and Carver 2012). Prior to the development of the MCS method simulations were used to examine definitive cases and statistical sampling was used to estimate the uncertainties in the simulations. MCS inverts this classic method, solving definitive problems by using a probabilistic analogy. This method of simulation is based on a computational algorithm that uses repeated random sampling to obtain its results and therefore requires (Singh and Arora 2015). According to , an MCS can be applied to complex systems in which variables can be adjusted simultaneously, and can be used to check for a synthetic effect by repeatedly sampling input values from their respective probability distributions. This is an essential approach for evaluating uncertainty when the uncertainty of the MCDA components can be characterized as a probability distribution or a confidence interval (Helton 2008;Janssen 2013;. Since the ANP technique determines criteria weightings in a network process on the basis of expert knowledge it is open to a variety of uncertainties. In order to overcome this particular problem, Rosenbloom (1997) proposed that the elements of a pairwise comparison matrix could be considered as random variables a i,j , with the proviso that it is in the distribution of 1/9 and 9, where a j,i = 1/a i,j and a i,i = 1. The pairwise comparison matrix is therefore symmetrically complementary. The numerical value of a random element a j,i will be related to a i,j , and it can therefore be assumed that {a i,j ji>j} is independent, and that the final scores S 1 , S 2 , …, S n will be stochastic. If S i > Sj, then alternative i is preferred to alternative j at a particular degree of error (a). In order to access information concerning the probability of element a i,j in the pairwise comparison matrix, it is assumed that all experts will produce the same evaluation of element a i,j . Every a i,j element will therefore be transferred into a discrete random variable. On the other hand, if there is only one decision maker, the compared value of each paired uncertainty will tend to be a continuous, random variable (Rosenbloom 1997; Hsu and Pan 2009;Feizizadeh et al. 2014b).
In this study, the ANP-MCS approach samples the vector of the input data in a random sequence in order to obtain a statistical sample of the vector of product variables. The ANP-MCS approach used to model the error propagation from the input parameters to the final results can be summarized by the following four steps: (1) A land subsidence inventory data-set is prepared as training data for the sensitivity analysis.
(2) Criteria weightings derived with the ANP method are used as reference weightings for the MCS (see Table 7). (3) The simulation is run N times with the number of simulations (N) recommended to be between 100 and 10,000), where N depends on three factors (a) the computational load, (b) the level of complexity of the model used, and (c) the accuracy required. (4) Spatial maps are produced of the results and the outputs analysed, with emphasis on the spatial distribution of the computed errors including the minimum rank (Figure 3(A)), maximum rank (Figure 3(B)), average rank (Figure 3(C)), and standard deviation rank (Figure 3 (D)). Figure 3(D) shows the spatial variation in uncertainty within the study area based on initial weightings derived from a traditional ANP. According to this Figure, in case of using the initial weightings of the ANP approach, the results in dark red spots will be more uncertain.

Variance-based global sensitivity analysis
Variance-based GSA is an approach that aims to measure the relative importance of different input variables (Guo and Du 2007). By adding error to the input data, the results of a model in which experts and/or researchers determine the initial values of the input variables can easily be affected, and this response can be measured by GSA. The values of criteria were determined on the basis of weightings derived for the model's input variables. It should be noted that such an ordering leads to an improved understanding of the indefinite parameters and to minifies the uncertainty interval of the response, which is necessary in order to obtain an acceptable uncertainty response range (Xu and Gertner 2008;Aven and Nøkland 2010;Hao et al. 2012). The results of the GSA determine those input variables that have the greatest effect on the results. (Ligmann-Zielinska and Jankowski 2012).
Variance-based methods have been applied to a wide range of problems; they are very helpful for sensitivity analyses and are appropriate methods to use with GSA Saisana et al. 2005;Feizizadeh et al. 2014b). Variance-based GSA can be used to determine which weightings have the greatest impact on the criteria ranking. Two sensitivity measures were produced by GSA, these being the first-order sensitivity index (S) and the total effect (ST) sensitivity index. According to Norton (2015), sensitivity analysis is often performed by looking at the output values affected by samples from a given distribution of parameter values. The GSA approach is a technique that assesses the interactive influence of a number of factors in a complex decision-making process. This method basically uses training data as samples for investigating the ways of addressing the experiment design problem of earning an acceptable compromise between sample size and the value of the data obtained from the outputs (Kleijnen 2005;Saltelli et al. 2008;Norton 2015). Since dealing with all outputs at a large number of points is both difficult and time-consuming, we define the sensitivities in terms of the expected values (probability-weighted averages) of some measure of variation, which makes it easier to understand and interpret outputs obtained using these expected values (Norton 2015). To this end, we fix some subset x 0 of a given input factor x and pay no attention to other factors that cause y to vary, where y denotes the output variable. Thus, y becomes y x 0 , say, with a mean of y x' , E [y x 0 ]. The fixing of x 0 is assumed not to bias the result. Consequently, The variance of y is defined as var(y) , E[ y À y ð Þ 2 ] because all influences on y can be easily divided between var(y x' ) and the mean of the variance remaining after x 0 has been fixed (Norton 2015): The formula (3) has been proved by Norton (2015). The variance of y x 0 from the decomposition of var(y) illustrates to what extent y x' can change the variability of y, or how strongly x 0 influences y. The variance ratio is therefore defined as: which can evaluate the sensitivity of y to x 0 . Since the ratio is independent of units, it can include any values between 0 (x 0 has no effect on y) and 1 (x 0 determines all the variability of y). The variancebased first-order sensitivity index of y to x j on its own can be obtained via the following equation: The variance-based first-order sensitivity index has been explained with the help of an example by Loucks et al. (2005). A broader measure of the significance of x j 's is the total sensitivity index S Tj of its effect on y, both on its own and through all interactions in which it is involved. If there are m input factors and x »j , then the m-1 parameters (all m parameters except x j ) have no interaction expect with x j ; S Tj is then the proportion of the var (y) that results from x»j alone has been accounted for, as follows (Norton 2015): For instance, in a model of m = 3 independent parameters, the 3 total sensitivity indices are defined as (Saisana et al. 2005;Feizizadeh et al. 2014b) S T1 ¼ S 1 þ S 12 þ S 13 þ S 123 and likewise For an input parameter x j , an important difference between S Tj and S j indicates a significant role for interactions with that parameter in y. Determining the interactions between given parameters enables a more comprehensive view of the model structure (Saisana et al. 2005;Feizizadeh et al. 2014b). A Sobol decomposition (Sobol 1993) was used to decompose the relationship between the resulting factor and the variance of the result. By using this approach allows sensitivities to be evaluated on the basis of the relative proportions of the resulting variance for each of the parameters, both separately and in combination. All parameters are considered to be independent, with values between 0 and 1. The relationship between the output y and a given parameter vector x is split into all components that have increased the number of their parameters (Norton 2015): Each component in ( Ã ) is given by etc. The variance of the result can be simplified through the orthogonality procedure: where i is the unit box of the range [0,1].
The sensitivity can be described as a variance ratio, as in (4); dividing (10) by V, , and so on (Norton 2015). The analysis was then continued to compute the significance of spatial bias in assessing the option rank order through the average shift in ranks (ASR), as follows (Saisana et al. 2005; Ligmann-Zielinska and Jankowski 2012): where a_rank ref is the rank of option A in the reference ranking (e.g. equal weight case), and a_rank is the last rank of option A. The ASR obtains the relative shift in the position of the entire set of options and quantifies it as the sum of absolute differences between the last option rank (a_rank) and the reference rank (a_rank ref ), divided by the total number of options (Ligmann-Zielinska and Jankowski 2012). In the first step of the analysis, weightings derived from the ANP are presented as the reference ranking column a in Table 8 and the results of the GSA presented in columns b, c, d, and e.

Results and validation
Two different approaches (ANP and SESUA) were used to produce LSSM maps. The weightings derived from each approach were applied through GIS spatial analysis and data aggregation models.
In the first analysis approach, the weightings calculated using the ANP technique were applied to produce a land subsidence susceptibility map. This approach was based on the common GIS-MCDA method for generating base maps, in order to be able to compare them with the results of the second approach. In the second approach, the LSSM map was produced by applying the results of the SESUA approach. Figure 4 shows the LSSM results for the two different approaches. In order to validate the results and measure the improvement in accuracy through a sensitivity analysis, a validation was carried out based on 23 GPS locations from the land subsidence inventory data-set. The validation step aimed to use statistically relevant algorithms to ensure the accuracy of the derived maps and determine their reliabilities (Sousa et al. 2004). The spatial effectiveness of both types of susceptibility map produced in this research was examined using an ROC curve, which is an effective method for characterizing the quality of deterministic and probabilistic prediction systems. This curve plots the false positive rate on the X axis and the true positive rate on the Y axis, thus illustrating the relationship between the two rates (Sezer et al. 2011). According to the theory behind the ROC curve, the area under curve (AUC) indicates the quality of a prediction system. This is determined by the system's ability to forecast the occurrence or non-occurrence of predefined 'events' (Negnevitsky 2002;Oh and Pardhan 2011). In order to obtain the relative ranks for each of the LSSM products, an index value of was calculated for each pixel of the resulting LSSM maps. All of these values were collected and organized in descending order. AUC values generally vary between 0.5 and 1.0, with the best results having an AUC value close to 1.0. A value close to 0.5 indicates that there is either no relationship or only a weak relationship that may exist purely by chance (Bui et al. 2012). The data that we used to draw up the ROC curve comprised 23 GPS coordinate sets for land subsidences that have occurred within the study area. The GPS coordinates for land subsidence locations were collected by the Ministry of Water Resources and Electric Power (MWREP) during the course of field work. This land subsidence inventory was used to validate the results with the ROC method. Figure 5 shows a comparison between the observed land subsidence locations and two LSSM maps of sub-areas within Marand County. ROC curves were plotted for both LSSM products and the AUC values calculated. In Figure 5, the rate of true-positives was plotted on the vertical axis, with the horizontal axis showing the rate of false-positives for both of the LSSM maps. The results obtained for the two LSSM products using the ROC method indicated an accuracy of 96% for the map derived using the SESUA approach, while the map derived using the ANP approach was assessed to have an accuracy of about 90%. The 6% improvement in accuracy clearly demonstrates that GSA is able to reduce the error in the ANP method and improve the accuracy of results. Figure 6 shows the results obtained using ROC and AUC methods to validate LSSM products against a land subsidence inventory data-set. Table 8 also presents the validation results for both of these approaches, together with the number of observed land subsidence events within the different ranges of land subsidence susceptibility from low (0%) to high (100%) and the area of each range. The high susceptibility range in the LSSM map produced using the ANP method contains only 17 of the 23 land subsidence events in the inventory data-set, with the remaining five land subsidence events falling in the moderate susceptibility range. Meanwhile the susceptibility map produced using the SESUA method yielded 19 high susceptibility locations and 4 with moderate susceptibilities (see Table 8). These figures indicate that the distribution of land subsidence events in the land subsidence inventory data-set fit better with the highly susceptible areas identified using the SESUA method than with those identified using the ANP method, thus demonstrating the improved accuracy of the ANP method when using the SESUA approach. The differences between the LSSM products obtained using both approaches can be seen in Figure 7, which presents enlargements from three different locations in each of the resulting maps.

Conclusion
This article has focussed on using the SESUA approach for sensitivity and uncertainty analysis of GIS-MCDA methods, by integrating the ANP, MCS, and GSA. Although GIS-MCDA-based methods are increasingly being used for spatial problems, their results are not necessarily reliable due to their uncertainties. Our study, which was part of an integrated GIS-based workflow, has increased the transparency of the results of GIS-MCDA methods by incorporating global sensitivity and uncertainty analyses. By adding these two analysis steps, the consequences of competing opinions obtained from a variety of experts and stakeholders are revealed and the influences of the chosen criteria become clear. While some authors have previously drawn attention to these problems in scientific publications (Ascough et al. 2008), the authors believe that only a workflow that is integrated within a GIS system will increase the accuracy of GIS-MCDA methods. This is achieved through the use of improved decision-making techniques based on error analysis methods, concepts, and assumptions. The implementation, monitoring and review, and possible correction of the selected methods (such as LSSM) will then become transparent, transferable and flexible. In this research we have been able to optimize the algorithms and achieve a degree of confidence for the resulting maps through use of the ROC method. As the results of this research have illustrated, further progress and optimization of the accuracy of GIS-based MCDA can be achieved through the use of spatially explicit methods. The precision of the weightings derived from pairwise comparison matrices obtained using the ANP method can be improved by incorporating GSA and MCS for sensitivity analysis.
In summary, these investigations have resulted in two main achievements, the first being the use of an SESUA approach to integrate GIS techniques with ANP techniques, and the second being the use of this approach to assess and minimize the sensitivity and uncertainty associated with LSSM, resulting in significant improvements in accuracy. Future research will focus on applying sensitivity and uncertainty analyses to data aggregation functions in a GIS environment, in order to assess and minimize the uncertainty of these methods. We will also optimize the ANP method itself by applying interval matrix methods. Results of this study are of considerable significance for the GIS community, and the authors are confident that this integrated methodology can be used in a variety of applications, and also for modelling purposes, within the GIS domain.