Presenting logistic regression-based landslide susceptibility results

(1)

Contents lists available atScienceDirect

Engineering Geology

journal homepage:www.elsevier.com/locate/enggeo

Presenting logistic regression-based landslide susceptibility results

Luigi Lombardo

a,b,⁎

, P. Martin Mai

b

a_{Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi}

Arabia

b_{Physical Sciences and Engineering (PSE) Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia}

A R T I C L E I N F O

Keywords:

Binary logistic regression Landslide susceptibility Standardized results

Least Absolute Shrinkage Selection Operator (LASSO)

A B S T R A C T

A new work-flow is proposed to unify the way the community shares Logistic Regression results for landslide susceptibility purposes. Although Logistic Regression models and methods have been widely used in geomor-phology for several decades, no standards for presenting results in a consistent way have been adopted; most papers report parameters with different units and interpretations, therefore limiting potential meta-analytic applications. Wefirst summarize the major differences in the geomorphological literature and then investigate each one proposing current best practices and few methodological developments. The latter is mainly re-presented by a widely used approach in statistics for simultaneous parameter estimation and variable selection in generalized linear models, namely the Least Absolute Shrinkage Selection Operator (LASSO). The North-east-ernmost sector of Sicily (Italy) is chosen as a straightforward example with well exposed debrisflows induced by extreme rainfall.

1. Introduction

The concept of landslide susceptibility dates back toBrabb (1984), being deﬁned as the likelihood of landslide occurring in a given area on the basis of local terrain conditions. For catchment, regional and larger scales (e.g.Shou and Lin, 2016) geomorphological mapping represents theﬁrst step to create inventories from which pseudo-quantitative or quantitative methods can be applied using a set of instability factors, in order to determine the landslide susceptibility for all mapping units (Hansen, 1984) partitioning the geographic space under study.

Pseudo-quantitative methods typically assign weights to predis-posing factor classes based on the frequency of landslides (e.g.Leoni et al., 2009). These weights are interpreted as the effect of each factor on landslide occurrence, while their sum over a given area roughly reflects landslide susceptibility. However, these ad-hoc empirical ap-proaches are unable to represent the likelihood (or probability) of landslide occurrence in a rigorous statistical sense, because they are not explicitly based on an underlying probability distribution. Conversely, with statistical methods (e.g.Ciurleo et al., 2017), the concept of sus-ceptibility is intrinsically linked to the probability of landslide occur-rence (in a give mapping unit), and therefore these quantitative methods offer more useful and detailed information by providing pre-cise results that are interpretable from a statistical perspective. Nu-merous statistical methods have been used in the landslide

susceptibility literature, ranging from Binary Logistic Regression (BLR) (e.g.,Rossi and Reichenbach, 2016), to non-parametric (e.g.,Guns and Vanacker, 2012) and Bayesian approaches (e.g., Lombardo et al., 2018a), each having their pros and cons. Among these methods, the most frequently applied is BLR, because of its wide applicability, ease of use, and software availability. Many geomorphological studies based on BLR have already been published, either on distinct geographic areas across the globe (e.g. Steger et al., 2016) or sometimes the same catchment (Cama et al., 2016), using a pixel-based discretization of space (Wang et al., 2013) or larger-scale slope units (Van Den Eeckhaut et al., 2009), focusing on landslide inventories triggered by heavy rainfall (Lin et al., 2017) or earthquakes (Parker et al., 2017), and they are often based on a diﬀerent but overlapping set of spatial predictors. Despite the great variability of modeling approaches and contexts, all these studies still share similarities that could potentially be exploited in reviews and subsequent meta-analyses, widening and reﬁning the perspective on slope responses to landslide triggers. However, there is currently no widely-adopted standards for presenting BLR results, hence drastically limiting the scope of such meta-analyses. For example, the comprehensive review ofBudimir et al. (2015)goes through 75 papers with applications across the globe, but because of the lack of consistency across papers for presenting results, this review is limited to general considerations, such as identifying and counting the predictors that are commonly used. Similarly,Pourghasemi and Rossi (2016)limit

https://doi.org/10.1016/j.enggeo.2018.07.019

Received 20 January 2018; Received in revised form 11 July 2018; Accepted 19 July 2018

⁎_{Corresponding author at: Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah University of Science and} Technology (KAUST), Thuwal, Saudi Arabia.

E-mail address:luigi.lombardo@kaust.edu.sa(L. Lombardo).

Available online 24 July 2018

(2)

their assessment to reporting the frequency of common covariates used in 220 papers. These qualitative considerations do not provide any quantitative information on the importance and eﬀect of each pre-dictor.

To facilitate future meta-analyses and quantitative comparisons across studies, we propose here a list of actions to be taken when modeling landslide susceptibility with BLR, in order to standardize the way BLR results are presented. Although the estimated susceptibility map and corresponding regression coefficients may depend on the choice of mapping units and general aspects that are dictated by the context, certain fundamental principles and diagnostic tools should be adopted more broadly in the geomorphology literature. Our proposed actions and work-flow, detailed inSection 2, concern six aspects of the whole statistical analysis procedure: i) Data description; ii) Variable selection; iii) Assessment of variable importance; iv) Landslide/no-Landslide proportion, implications and probability cutoff choice; v) Predictive performance; vi) Single and multi-fold modeling procedures. Our proposed framework is illustrated throughSection 4, using as a practical example the Multiple Occurrence of Regional Landslide Event (MORLE,Crozier, 2005) that took place in North-East Sicily, Italy, on October 1, 2009. This choice is justified due to the almost unique landslide class that was triggered, namely debrisflows. This makes the implementation and interpretation easier also due to the large number of papers published on this specific MORLE.

2. Overview of various practices in landslide susceptibility

Data description: Thefirst discrepancies among articles can be found in the way the coefficients of the linear covariates are published, being standardized only in few cases (e.g. Raja et al., 2017) or simply re-ported in the majority of contributions (e.g.Erener et al., 2016), if re-ported at all. The very fact that covariates are generally kept in their original units hinders the comparison of their impact with respect to other covariates in the same study area and the comparison of the same covariate in different geographic contexts. Enabling the latter could open interesting discussions and generalized interpretation on landslide activations. Variable selection: Variable selection procedures should be always included in any work-flow for landslide susceptibility to remove non-informative predictors, avoid model over-complexity and variable-interaction issues and deal with multicollinearity (Marchesini et al., 2014). This is another case where strong differences arise among re-search groups. In some cases, no variable selection takes place (e.g., Feng et al., 2016) whereas in the opposite situation, the community almost unanimously performs a Stepwise Selection (e.g., Eeckhaut et al., 2010). However, the statistical community has highlighted nu-merous deficiencies for such method since 1990's (Copas and Long, 1991;Derksen and Keselman, 1992;Harrell Jr, 2015) and few cases are currently available where more effective variable selection procedures are tested in landslide susceptibility studies (e.g.Castro Camilo et al., 2017). Assessment of variable importance: Almost no agreement can be found within the community when assessing the strength of each cov-ariate with respect to landslide occurrences. Few cases take into con-sideration the p-values for significance (Chen and Wang, 2007), other articles” multiply the maximum value of the variable in the dataset with its regression coefficient” to obtain a metric defined as Measure of Parameter Importance (MPI, Guns and Vanacker, 2012) and others consider the frequency and rank of selection in a stepwise framework (Cama et al., 2017). Landslide/noLandslide proportion, implications and probability cutoff choice: Very few articles explore the effects of the ba-lanced or non-baba-lanced sampling scheme used to construct the di-chotomous dataset to be fitted (Heckmann et al., 2014; Conoscenti et al., 2016). Even more rarely, we assess the implications that such choice produces onto the susceptibility itself (Petschko et al., 2014) and onto the cutoff value that separates probabilities into stable and un-stable predicted conditions (Frattini et al., 2010). Predictive perfor-mance: A better agreement must also be found on model performance

descriptions across papers. Most of the community uses Receiver Op-erating Characteristic (ROC) curves and their Area Under the Curve (AUC); however, not always a better insight is provided through con-fusion matrices and their implications (e.g.Rossi et al., 2010). In this regard, numerical quantification of performances are clearly funda-mental although reporting ROCs and AUCs only, does not yield the geological information of the model hits (True Positives and True Ne-gatives) and misses (False Positives and False NeNe-gatives). Single and multi-fold modeling procedures: This overall framework is further di-versified considering that some papers merely regress the landslide data once (e.g.Lee, 2014) and publish the corresponding results, while other authors implement more rigorous multi-fold procedures to describe estimated model parameters through their mean and variance based on a number of independent bootstrap replicates (Brenning et al., 2015). This is the case for the K-fold cross-validation procedures described by Petschko et al. (2014).Petschko et al. (2014)also mention something of great relevance which we will further investigate in this contribution. They highlight the differences between what we will define as “true” landslide occurrence probabilities in contrast to the concept of landslide “susceptibility”. The way most of the community computes landslide susceptibility maps, by selecting a priori a balanced number of landslide presences and absences, provides biased estimates of the“true” pixel-based probabilities. We here clarify this concept through synthetic numerical simulations (see Supplementary Materials), and show that, although susceptibility maps look similar in both cases and conclusions about slope stability are essentially the same,“true” probabilities are not a linear function of“susceptitiblity” values, and using the former improve result interpretation over the latter.

3. Study area and landslide inventory

The adjacent catchments of Scaletta and Itala are located on the eastern side of the Peloritani Belt (Fig. 1b) and included in the center of the 1st October 2009 storm that caused 37 victims and 0.5 M€ worth of structural damage. At the time of the disaster the area including the two catchments was exposed to 75 mm and 23 mm of rain two and one weeks before the main storm which discharged 250 mm in 8 h (Bout et al., 2018). As a result, the saturated weathered cover draping over a parent bedrock of medium to high grade metamorphic rocks (Fig. 1c) destabilized and gave rise to hundreds of shallowflow-like landslides. These debris flows have been mapped on the basis of Google Earth (2006–2010) and orthophotos (six days after) provided by the National Civil Protection (Cama et al., 2015). Pre- and post- event debrisflows have been intersected to isolate the component triggered by the storm in 2009.

The total inventory (Fig. 1d) for the two catchments accounts for 1098 debris flows. A histogram and related statistics of the corre-sponding maximum distance from the crown to the runout zone is shown in Fig. 2 where it emerges that most debris flows traveled for > 40 m, indicating high water content andfluidity of the moving mass.

4. Statistical modeling

4.1. Dataset creation

Several representations of the geographic space have been proposed throughout the years. These are deﬁned as mapping units (Carrara et al., 1995) and represent the spatial object partitioning a given study area. Among various mapping units, we adopt a 2 m-side squared lattice coinciding with the resolution of the Digital Elevation Model used for the analyses. The typical approach for building datasets also includes designing the unstable regions for each landslide triggering location. Diﬀerences in the literature can be found in groups representing land-slide presences as the centroid of the mapped polygons (Hussin et al., 2016) or as the highest pixel along a landslide crown (Landslide

(3)

Identiﬁcation Point or LIP,Lombardo et al., 2014). In addition to this, authors have suggested to widen the instability regions around the aforementioned locations. This is the case for“seed” or “diagnostic” areas proposed by Süzen and Doyuran (2004)and Rotigliano et al. (2011), respectively. The reason behind the two aforementioned ap-proaches is to model landslide instabilities as an areal rather than a point process. However, in this contribution we avoid seed and diag-nostic areas for they may give rise to highly dependent datasets which in turn will conﬂict with the independence assumption of the likelihood function described inSection 4.3. As a result the number (1098) and locations of unstable pixels here coincides with the LIPs. In subsequent statistical modeling and analysis, a balanced sample of landslide ab-sences will be randomly extracted among the remaining stable pixels.

4.2. Covariates

Fourteen covariates are used to explain the spatial distribution of landslides, based on the dataset obtained as explained inSection 4.1.

Two covariates are categorical and correspond to the i) Outcropping Lithology (Lentini et al., 2000)1and ii) Land Use (Corine 20062). The remaining twelve covariates are continuous, nine of which are obtained from a pre-event 2 m h-resolution and 0.25 m z-resolution LIDAR Di-gital Elevation Model (ARTA 2008, Cama et al., 2015). These nine covariates are: iii) Elevation, iv) Eastness and v) Northness (Lombardo et al., 2018b), vi) Plan and vii) Proﬁle Curvatures (Heerdegen and Beran, 1982), viii) Relative Slope Position (Böhner and Selige, 2006), ix) Slope (Zevenbergen and Thorne, 1987), x) Stream Power Index (Moore et al., 1991), xi) Topographic Wetness Index (Beven and Kirkby, 1979). The 15 m resolution xii) Normalized Diﬀerence Vegetation Index (NDVI) (Rouse Jr et al., 1974) is calculated on the basis of an ASTER scene acquired on July 2007. Distances to xiii) Fault Lines and to xiv) Channel Network are computed as Euclidean distances from every pixel Fig. 1. Italy and focus on study area (a); Itala and Scaletta catchments (b); Outcropping Lithologies and subcatchments (c); Landslide inventory, 1st October 2009.

1_{http://www.isprambiente.gov.it/Media/carg/601_MESSINA_REGGIO/} Foglio.html

(4)

of the 2 m lattice to the nearest tectonic alignment and channel. All the covariates are finally resampled to the 2 m resolution. Their ID, ac-ronym and classes for categorical variables are shown inTable 1, they will be alternatively used infigures through the manuscript. Ultimately, all the continuous covariates are rescaled by substracting their mean value and dividing by their standard deviation. This last step is crucial to be able to objectively assess the relative importance of each covariate by comparing their corresponding estimated regression coefficients. It might, however, be easier to interpret regression coefficients on the

covariates' original scale. A simple linear transformation can be used to get the regression coeﬃcients on the original scale from the ones ob-tained on the normalized scale; this is further detailed in Appendix A.

4.3. Parameter estimation and variable selection via the LASSO

Details on Binary Logistic Regression models are provided in AppendixA whereas this section will focus on how parameters can be estimated, and how to perform variable selection using the LASSO (Tibshirani, 1996).

Suppose that n independent observations y1,…, yn∈{0, 1} from the

logistic model(5)are available. In our context, these observations re-present both landslide presences and absences. For each i = 1,…, n, denote byxi⋆= (x1i⋆,…,xpi⋆) the vector of normalized covariates

as-sociated to the observation yi, and byπ(xi⋆) the probability that the ith

observation equals one. Regression parameters (on the normalized scale) may be estimated by maximizing the log-likelihood function, which may be written as

∑

… = + − ⋆ ⋆ ⋆ = ⋆ = ⋆

(

β β β

)

π x π x ℓ , , , _p log{ ( )} log{1 ( )}, i y i i y i 0 1 :i 1 :i 0 (1)

whereπ(xi⋆) is linked to the regression parameters through(5). E

ﬃ-cient algorithms to maximize [1] have been developed, and are freely available in the statistical software R.

When the number of covariates p is large compared to the sample size n, when many of these covariates are not signiﬁcant, or when they are highly multicollinear, the performance of estimators might be highly aﬀected and strategies to achieve variable selection need to be found. One convenient and powerful variable selection technique, which can be viewed as an alternative to stepwise selection methods, performing simultaneous estimation and variable selection is the Least Absolute Shrinkage Selection Operator (LASSO). Essentially, the LASSO maximizes a penalized counterpart of [1], which may be expressed as

∑

… = + − − ⋆ ⋆ ⋆ = ⋆ = ⋆ = ⋆

(

β β β

)

π x π x λ β ℓλ , , , p log{ ( )} log{1 ( )} , i y i i y i j p j 0 1 :i 1 :i 0 1 (2) whereλ > 0 is the penalty parameter, which is usually chosen by cross-validation in order to optimize prediction accuracy. Whenλ = 0, (2)boils down to(1), and whenλ → ∞, all estimated parameters (ex-cept the inter(ex-cept β0⋆) are set to zero, i.e., = ⋯= =

⋆ ⋆

β₁ β_p 0. In be-tween, it can be shown, thanks the geometry imposed by the penalty term in(2), that some parametersβj⋆may be set to zero while others

might be non-zero; by setting some parameters to zero, this auto-matically achieves variable selection.

4.4. Implementation

We write the code using the statistical software R (seehttps://www. r-project.org/) and use the glmnet package written byFriedman et al. (2009)to implement the LASSO. The adoption of a 2 m-grid mapping unit partitions the study area into 3.8 · 106pixels. To maintain a ba-lanced sampling scheme throughout the analyses we i) extract from the dataset all the cases with a status assigned as 1 (landslide presences) and ii) merge an equal number of cases with status assigned as 0 (landslide absences) randomly selected from the remaining data. We then randomly split this smaller dataset into two balanced subsets, each containing 75% and 25% of the cases, respectively. Thefirst subset (with 75% of cases) is used forfitting whereas the second (with 25% of cases) is used for validation. This operation is then repeated 500 times to producefitted models. For each replicate, the value of the penalty parameterλ in (2)is re-estimated by 10-fold cross-validation, and a different number of variables is selected by the LASSO. For each of the 500 modelsfitted, the regression coefficients βjare stored generating a Fig. 2. Distribution of runout distance calculated for each debrisflow from the

crown to the deposition area where possible, otherwise to the intersection with the river network.

Table 1

Continuous and categorical predictor list. Throughout the manuscript, IDs are used inﬁgures while acronyms and classes are used in the text.

ID Acronym Class 1 Dist2Faults – 2 Dist2Channel – 3 Eastness – 4 Elevation – 5 NDVI – 6 Northness – 7 Plan_Cur – 8 Prof_Cur – 9 RSP – 10 Slope – 11 SPI – 12 TWI –

13 Use111 Continuous Urban Fabric 14 Use323 Sclerophyllous Vegetation

15 Use314 Partially Wooded Land or Degradated Forest 16 Use221 Vineyards

17 Use231 Pastures

18 Use322 Moors and Heathland 19 Use223 Olive Groves 20 Use331 Beaches, dunes, sands 21 Use321 Natural Grassland 22 PMAa Paragneiss to Micaschists 23 FDNa Muscovite Marble 24 PMPa Pegmatite/Aplite

25 MLEa Crenulated Paragneiss/Micashists 26 FDNb Phyllites/Meta-Arenites 27 MLEc Double Micas Marbles 28 VEP Metarenite/Metapelite 29 PMAd Silicate Marble

30 bn4 Terraced Alluvial Deposit (4th) 31 bn3 Terraced Alluvial Deposit (3rd) 32 gn2 Terraced Marine Deposit (old) 33 b2 Eluvium/Colluvium 34 bb Alluvial Plain Deposits 35 gn3 Terraced Marine Deposit (recent)

(5)

matrix B of dimension 500 by 14, where columns denote covariates and rows denote replicates. To assess variable importance, we then compute the proportion of times each covariate was selected, i.e., when the corresponding regression coeﬃcient βj⋆ was not set to zero by the

LASSO. For each of the 500 replicates, we also compute the proportion of True/False Positives/Negatives, the Receiver Operating Character-istic (ROC) curve and Area Under the Curve (AUCs), the error rate (mean(FP|FN)), the ﬁtted probabilities and the response values (sus-ceptibility V S covariates) for each of 3.8 · 106_pixels.

In addition we estimate the best probability cutoff between stable and unstable conditions (Castro Camilo et al., 2017). This is done by calculating the Accuracy (TP + TN) divided by the Total Population (TP + TN + FP + FN) at various probability cutoffs. This operation ensures that the cutoff choice, meant to predict landslide presences (ones) and absences (zeros) from estimated probabilities, is quantita-tively determined rather than subjecquantita-tively adopted a priori. The idea behind this test comes from the considerations inFrattini et al. (2010) andHeckmann et al. (2014)where authors stress the fact that prob-ability cutoffs are data-dependent. Here we select the cutoff that max-imizes on average the Accuracy across 500 replicates.

5. Results

Fig. 3a demonstrate how we propose to choose the best probability cutoff by selecting the maximum accuracy value along the dashed curve expressing the average across 500 random replicates. The community usually adopts a 0.5 cutoff which is also reflected in the present case confirming its correctness for balanced sampling strategies. Fig. 3b shows 500 ROC curves and their overimposed mean corresponding to an AUC of 0.84 and a standard deviation of 0.01. As the ROC/AUC provides just a lumped measure of model performances, this informa-tion is complemented in Fig. 4a with a scatterplot representing the combination of TP/Observed Positives and TN/Observed Negatives for each of 500 models.

Even in this case a limited variability is shown with standard de-viations of 0.80 and 0.86 along the abscissa and ordinate, respectively. The mean proportion of TN and TP is 79.70% and 75.02%, respectively. Considering these central tendencies and their associated variabilities it emerges a better prediction toward stable conditions.

Error rates inFig. 4b conﬁrm the overall good performances with a median value of 0.228 and an inter-quartile distance of 0.008.

Once a strong predictive power and a low variance among replicates are proven, the geological interpretability of the results can be assessed. We do this in the domain of the covariates and final probabilities. Specifically, in this section we combine the results with a commentary on the main novelty we propose for presenting the covariate effects in a way that the community can compare even for different sites and events.

For covariates, Fig. 5 summarizes the procedure described in

Section 4.4. The percentage of variable selections is shown in the top panel in a colour coded bar plot with assigned red and blue for positive and negative betas, respectively. The sign is assessed at the mean of each of the 500 coefficients. The lower panel shows the absolute value of the standardized beta coefficients. As the beta are in the same unitless scale it is easy to recognize dominant covariates. The inter-pretation of the standardized betas still indicates that the probability increases or decreases by a factor equal to the exponential of the beta value. Whether that would increase or not the probability is lead by the sign colorcoded above. We also include the percent contribution to the model (samefigure, second axis to the right). This is done by normal-izing the rescaled betas with respect to the maximum beta value.

Fig. 6shows the results of Jackknife tests (Lombardo et al., 2016b). Jackknife procedures generate models with one variable only to be associated with the corresponding models built with all variables but one. Each model undergoes the same K-fold resampling scheme of the LASSO-penalized models and their AUCs are compared with a reference full model made by using all available variables. This allows one to i) highlight the univariate predictive power of single covariates; ii) stress the same idea but in a multivariate environment where the importance of a given covariate is shown in terms of performance drops. This re-presentation here stresses the role of the Slope. This would have ap-peared to be less relevant than the Elevation only looking at the betas in Fig. 5. However, both Jackknife tests indicate the Slope to be able to describe the landslide scenario the most. The Slope-only-model has an AUC equal to 0.64 whereas the all-but-Slope model reports the greatest drop with aΔAUC with respect to the full model equal to 0.01.

Another way to assess the effect of the covariates consists of studying the bi-variate relations between the spatial distribution of each covariate against the spatial distribution of thefitted probabilities (Tziritis and Lombardo, 2016;Schillaci et al., 2017).Fig. 7summarizes the responses of the covariate mostly selected by LASSO. Due to the different nature of covariates, for continuous cases we used a 2D kernel density estimator whereas for categorical variables we created an his-togram extracting the probabilities at locations with categorical values equal to 1. The interpretation of these plots indicates variations in the debris flow susceptibility in a multivariate domain. This means that dominant variables should reflect their effect while non significant covariates may not have enough influence to produce a change over the fitted probabilities.

Few demonstrations of the concept above are shown for Distance to Faults and Slope where with opposite effects they clearly contribute to modify the susceptibility despite the co-existing effects of the other covariates on the susceptibility. As for the single dummies, lithology bb (Alluvial Plain Deposit) has a negative beta sign and this is reflected in most of the probabilities centered around 0.1. A non-influential ex-ample can be seen in land use 221 (Non Irrigate Arable Land) for its beta sign is negative but the susceptibility values have an approximate normal distribution centered at 0.4 and a standard deviation 0.19.

(6)

The susceptibility map is calculated as the meanfitted probability per pixel across 500 replicates andFig. 8presents it (a) along with the corrected probabilities (b) and their respective standard deviations (c and d) measured in two intervals. Susceptibilities do not strictly cor-respond to probabilities in a statistical sense. This is due to the effect of the balanced sample and a detailed discussion can be found in (Petschko et al., 2014). We here report both the distorted probabilities corresponding to the susceptibility and their corrected equivalent fol-lowing the procedure described by Petschko and co-authors. In addi-tion, we extend the discussion from the aforementioned authors by comparing the effects of balanced and full samples in a simulation study presented in the Supplementary Materials. InFig. 8we reclassify the corrected probabilities into 5 classes corresponding to 2.5, 25, 50, 75 and 97.5 percentiles. These values are calculated from the true prob-ability distribution extracted at Landslide Identification Point locations in accordance to the suggestions inPetschko et al. (2014). The only difference is that Petschko and co-authors use only three quantiles (low, medium and high probabilities) rather than five (low, low-medium, medium, medium-high and high). Despite the patterns in the mean

predictive maps are similar, those in the standard deviation maps are not. This is due to the fact that the variance of the corresponding odds computed with the full model or V ar(odds) are not linearly linked to the variance of the odds obtained with a balanced sample which are mul-tiplied by a constant c. Hence, the latter becomes V ar(c · odds) = c2_{· V}

ar(odds).

Susceptibility and 95% confidence interval maps provide visual support for assessing the relative likelihood for debrisflow activations. However, a better insight on the relation between susceptibility and its stability across replicates is achieved through error plots (Guzzetti et al., 2006;Lombardo et al., 2016a).Fig. 9shows the bi-variate dis-tribution of the two aforementioned maps. We opt for a 2D kernel density estimator for a better visualization. Most of the variability is confined below 0.05 attesting for coherent replicates despite the 500 resampling iterations. Three peculiar trends emerge from the density plot. These are also shown in (Lombardo et al., 2015) and are due to categorical classes highlighted in the map in three areas, two of which toward the East and close to the coastline while the remaining one is shown to the West side of the main stream (seeFig. 8).

Fig. 4. TP proportion VS TN proportion (a). Red lines correspond to mean values whereas blue and green lines represent one and two standard deviation in-tervals, respectively. Error Rates (b). Both cross-va-lidation plots account for 500 random replicates. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web ver-sion of this article.)

Fig. 5. Frequency of LASSO selection across 500 random replicates. Red indicates a positive value for the mean beta while blue represents a mean negative correlation with respect to landslide presences. Variable labels are also valid for lower plots (1st row) and their IDs are shown in the x-axis. Beta coefficients calculated for rescaled continuous covariates and they are not modified for categorical classes. Orange line represents non-significance (2nd row). (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)

(7)

6. Discussions

What we can infer from the present example is that with respect to the average debris ﬂow susceptibility model Slope (β = 0.53),SL

Elevation (β_EL=−0.71), Distance to Faults (βD F2 =−0.44), Northness

(β_NN= 0.40) dominate the susceptibility distribution from the four diﬀerent perspectives described above among continuous covariates. This can be respectively interpreted as i) a clear eﬀect of the slope over shallow landslides; ii) topographic control over the rainfall pattern, which actually got dissipated mid-way between the outlet and the upper catchment divide; iii) fragmentation of the material in proximity to tectonic features, consequent release of material from common ero-sion processes and remotion of potential sources for shallow landslide; iv) solar radiation control over evapotranspiration processes and con-sequent retainment of soil moisture in North-facing slopes.

With respect to categorical cases, lithology bb or Recent Alluvial Deposits (β = −2.42), MLEc or Double Micas Marbles (βbb MLEc=−2.14)

negatively impact the likelihood whereas FDNb or Phyllites to Meta-Arenites (βFDNb= 0.47) increases the probability. We focused on these

lithotypes as they are frequently included, present a significant effect and widely outcrop in the area bb = 3%, MLEc = 36%, FDNb = 33%. From an interpretative standpoint their effect can be justified due to geotechnical and geomoprhological controls over instability. Phyllites in the area are strongly weathered and their foliated nature could give rise to preferential slide zones from micro to macro scales where the main triggering movement could take place. As regards the other al-luvial deposits and marbles, thefirst one typically occupies flat stable areas whereas the second is a competent rock by nature and thus it is not easy to mobilize it as part of a shallow landslide process.

Land uses appear to have a limited inﬂuence with the exception of Use231 or Pastures (βUse231= 0.50) and Use321 or Natural Grassland

(βUse321= 0.93). These both positively contribute to slope instabilities

and respectively outcrop for 12% and 7% of the area. They both in-dicate conditions where no dense vegetation can shield the soil from direct rainfall impact.

The interpretation can be further improved looking at response plots (Fig. 7). The eﬀects of Elevation, Slope, Distance to Faults clearly

emerge on the probability despite its dependency from the multivariate approach. The best example is that of the Elevation where the sus-ceptibility increases up to its maximum coinciding with an height of approximately 500 m.a.s.l. These values are lead by other covariates as they cannot be justified from the Elevation (negative beta). However, from the 400 m nick-point to the right tail of the distribution, the sus-ceptibility sharply decreases due to the negative effect of the Elevation. Mean susceptibility and error maps (Fig. 8) together with error plots (Fig. 9) are known to belong to level 4 assessment (Guzzetti et al., 2006). We extend this assessment to account for the probability cor-rection described byPetschko et al. (2014). This topic is not trivial as the distortion of the balanced sample only allows for a visual inter-pretation in terms of relative likelihood whereas a corrected probability is interpretable as the probability of a landslide under similar triggering conditions. Furthermore, the probability distortion affects the variance computed for each pixel in a k-fold cross-validation scheme, resulting in a strong overestimation.

From a geographic perspective, in the present study we modeled part of the area investigated byCama et al. (2015)but including the adjacent catchment to the North (Scaletta Zanclea). Despite the dif-ferent implementation, mean susceptible areas well coincide with the self-validated model 2009 by Cama and co-authors whilst the disper-sion they encountered is much higher. This can be due to two different reasons. Thefirst one being a different number of replicates as we run 500 random extractions instead of 10. The second reason could be the different variable selection procedure. Stepwise selection is not as strong as the LASSO due to the fact that the AIC criterion tends to in-clude more variables while LASSO operates on the domain of the beta coefficients setting every negligible contribution to a 0 coefficient.

The advantages of implementing LASSO can also be exploited to improve other susceptibility models. In fact, LASSO is not limited to statistics and has recently been successfully combined to data mining methods such as Neural network, Classiﬁcation And Regression Tree and Support Vector Machine (SVM) (e.g.Goo et al., 2016). Because engineering geologists and geomorphologists are increasingly using data mining to generate predictive models for landslides, LASSO could represent a valid tool too boost performances and reduce over-complexity at the same time.

Fig. 6. Jaskknife tests for continuous and categorical covariates. Only-one-variable models are shown in the 1st row where red lines correspond to a random prediction performance or 0.5 AUC. All-but-one-variable models are shown in the 2nd row where green lines are the average prediction of the full model. Blue lines are the average prediction of 500 LASSO penalized replicates. (For interpretation of the references to colour in thisﬁgure legend, the reader is referred to the web version of this article.)

(8)

(9)

Ultimately, the work-ﬂow we propose can be used in any context where landslides are modeled over space, from catchment to territory-wide scales and from susceptibility to hazard models. Even for landslide hazard studies (Ko and Lo, 2018), where the spatial component is in-itially disentangled from the temporal one during the analyses, LASSO can contribute to simplify the model whereas the analyses on the cov-ariate role and the validation procedures we present can boost the in-terpretation of the results from a geological and engineering

perspective.

7. Conclusion

The community has long invested eﬀorts in developing landslide susceptibility models. However, no clear standards are still in place with respect to some key part of the analyses. In terms of methodolo-gical development we present the LASSO penalization instead of its common Stepwise counterpart. Stepwise selection is subjected to a wide spectrum of disadvantages and technical problems. This is common knowledge within the statistical community and some the most com-prehensive description of these issues can be found inHarrell Jr (2015) orMiller (2002).

Taking aside the technical considerations, we propose a compre-hensive way to present any results for landslide susceptibility studies. Part of the suggestions were already summarized by Guzzetti et al. (2006) and subsequently adopted by other authors (e.g. Lombardo et al., 2015). However, we here propose an update. Thefirst step should always control the correctness of the probability cutoff. This cutoff is subjected to variation due to the sampling scheme and data as already reported byFrattini et al. (2010) andHeckmann et al. (2014). The maximization of the accuracy against a varying cutoff ensures an op-timal and quantitative cutoff selection irrespective to the data.

The second consideration relates to the way the community oper-ates in a cross-validation framework. True and False Positive propor-tions are not always included and this even more common when con-sidering cross-validated results. However, the information conveyed by Fig. 7. Bivariate distribution between most selected continuous covariates and the corresponding susceptibility value for each mapping unit. Categorical covariates are presented by summarizing the univariate distribution of the speciﬁc class if present, for each mapping unit. Red lines separate predicted stable from unstable conditions. (For interpretation of the references to colour in thisﬁgure legend, the reader is referred to the web version of this article.)

Fig. 8. Mean (a) and 95% conﬁdence interval (c) maps calculated across 500 random susceptibility replicates. Corrected and reclassiﬁed probabilities (b) and their standard deviation (d) are indicated by the asterisk. Dashed lines represent contours lines expressing the original density of landslides per meter squared whereas the over-imposed numbers are the actual densities.

Fig. 9. Error plot representing the bivariate distribution between mean sus-ceptibility and its 95% conﬁdence interval.

(10)

confusion matrices complements the compressed information from ROC and AUC which only describe performances from a numerical per-spective.

The third and main message of this contribution relates to the way the community shares logistic regressed coefficients and/or covariate effects on models. As things currently are, there is very little chance to compare different coefficients from a paper to another. However, this is done in other scientific branches supporting meta-analytic studies. Meta-analysis within the geological community could improve our understanding of environmental influences over different landslide types and triggers. Few authors already present coefficients after stan-dardizing the original covariates and we stress this to be the ideal choice to be additionally expressed in absolute values. We also compute the percentage contribution, count the frequency of inclusion, perform Jackknife tests and response plots to complement the covariate influ-ence assessment. Doing this could open new interpretation on covariate effects. A very simple but effective comparison across articles could be achieved by normalizing the absolute betas as shown in the 2nd row of Fig. 5to get an objective influence expressed in percentage. Further-more, Jackknife tests are convenient ways to express covariate influ-ences into AUC values, these being common metrics that every member

of the community knows and recognizes.

The overall work-ﬂow is intended to provide a guideline for pre-senting susceptibility results and we believe it could improve the quality of several papers and more importantly the cross reference/ comparison among contributions, ultimately allowing for meta-analytic applications or reviews even in landslide susceptibility practices.

Ultimately we investigate the eﬀect of a balanced sample over the probabilities extending the work of Petschko et al. (2014). Supple-mentary Materials include simulations and subsequent comparison between a susceptibilities and corrected probabilities.

Acknowledgement

The authors would like to thank Dr. Daniela Castro Camilo as the code used throughout the analyses is a slight modification of the LUDARA code included inCastro Camilo et al. (2017). Part of the sa-tellite images used to generate the landslide inventory were obtained thanks to the European Space Agency Project (ID: 14151) titled: A re-mote sensing based approach for storm triggered debrisflow hazard modeling: application in Mediterranean and tropical Pacific areas. Principal Investigator: Dr. Luigi Lombardo.

Appendix A. Statistical modeling using binary logistic regression

In classical regression analysis, we seek to explain the variability of a response variable Y in terms of p covariates x1,…, xp, such as the ones

presented inSection 4.2. When such covariates are assumed to be continuous for simplicity, the ordinary Gaussian linear model may be expressed as

= + + ⋯+ +

Y β₀ β x_{1 1} β x_{p p} ε, ₍₃₎

where β0, …, βp are regression coeﬃcients measuring the eﬀect of each covariate, ε ∼ N(0, σ2) is the random Gaussian component, and

β0+β1x1+⋯ + βpxpis the systematic component, which determines the expected value of the response Y. This framework can be extended in

many ways, e.g., by considering categorical covariates, non-linear effects, or non-Gaussian responses. Generalized linear models (GLMs), which have been studied in depth and widely applied in a variety of scientific fields, can handle responses whose distribution belongs to the linear exponential family, including the normal, binomial, bernoulli, Poisson, or gamma distributions, among others.

When Y has two possible outcomes 0 and 1 (representing, in our context, landslide presence/absence in a given pixel), a generalized linear model may be formulated as follows: denoting byπ the probability that Y takes the value 1, one sets

= + + ⋯+

η π( ) β₀ β x_{1 1} β x_{p p}, ₍₄₎

whereη is the link function relating the linear predictor on the right-hand side of(4)toπ, i.e., the mean of Y. In other words, Y follows a Bernoulli distribution with parameterπ = η−1₍_β

0+β1x1+⋯ + βpxp). Several link functions may be chosen. A common choice with nice theoretical

prop-erties is the logit functionη(π) = log {π/(1 − π)} with inverse link η−1(x) = exp (x)/{1 + exp (x)}, leading to Binary Logistic Regression (BLR). It can be veriﬁed that the odds, i.e., the ratio Pr(Y = 1)/Pr(Y = 0) = π/(1 − π), increase by a factor exp(βj) for a unit increase of the covariate xj, j = 1,

…, p. In particular, in the landslide context, when βj> 0, the covariate xjhas positive eﬀect on the landslide occurrence probability, while it has a

negative eﬀect when βj < 0.

Notice that, using the centered and rescaled variables xj⋆= (xj− mj)/sj, j = 1,…, p, with the mjs and sjs respectively denoting the mean and

standard deviation of each covariate xj, Eq.(4)may be rewritten as

= + + +⋯+ + = + +⋯+ + +⋯+ = + +⋯+ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ η π β β s x m β s x m β β m β m β s x β s x β β x β x ( ) ( ) ( ) ( ) , , p p p p p p p p p p p 0 1 1 1 1 0 1 1 1 1 1 0 1 1 ₍₅₎

whereβ0⋆=β0+β1m1+⋯ + βpmpandβj⋆=βjsj, j = 1,…, p. Hence, the coeﬃcients βjon the covariates' original scale can be recovered from the

coeﬃcients βj⋆on the normalized scale using the relations in(5), and vice versa. Coeﬃcients βjon the original scale are usually easier to interpret,

while coeﬃcients βj⋆on the normalized scale allow to compare the relative importance of each covariate on the response.

Appendix B. Supplementary data

Supplementary data to this article can be found online athttps://doi.org/10.1016/j.enggeo.2018.07.019.

References

Beven, K., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant. Hydrol. Sci. J. 24 (1), 43–69.

Böhner, J., Selige, T., 2006. Spatial prediction of soil attributes using terrain analysis and climate regionalisation. Gottinger Geographische Abhandlungen 115, 13–28.

Bout, B., Lombardo, L., van Westen, C., Jetten, V., 2018. Integration of two-phase solid fluid equations in a catchment model for flashfloods, debris flows and shallow slope failures. Environ. Model Softw. 105, 1–16.

Brabb, E.E., 1984. Innovative approaches to landslide hazard and risk mapping. In: 4th International Symposium on Landslides, 16–21 September, Toronto, Canada.

Brenning, A., Schwinn, M., Ruiz-Páez, A., Muenchow, J., 2015. Landslide susceptibility near highways is increased by 1 order of magnitude in the Andes of southern Ecuador, Loja province. Nat. Hazards Earth Syst. Sci. 15 (1), 45–57.

(11)

Budimir, M., Atkinson, P., Lewis, H., 2015. A systematic review of landslide probability mapping using logistic regression. Landslides 12 (3), 419–436.

Cama, M., Lombardo, L., Conoscenti, C., Agnesi, V., Rotigliano, E., 2015. Predicting storm-triggered debrisﬂow events: application to the 2009 Ionian Peloritan disaster (Sicily, Italy). Nat. Hazards Earth Syst. Sci. 15 (8), 1785–1806.

Cama, M., Conoscenti, C., Lombardo, L., Rotigliano, E., 2016. Exploring relationships between grid cell size and accuracy for debris-ﬂow susceptibility models: a test in the Giampilieri catchment (Sicily, Italy). Environmental Earth Sciences 75 (3), 1–21.

Cama, M., Lombardo, L., Conoscenti, C., Rotigliano, E., 2017. Improving transferability strategies for debrisﬂow susceptibility assessment: application to the Saponara and Itala catchments (Messina, Italy). Geomorphology 288, 52–65.

Carrara, A., Cardinali, M., Guzzetti, F., Reichenbach, P., 1995. GIS technology in mapping landslide hazard. In: Geographical Information Systems in Assessing Natural Hazards. Springer, pp. 135–175.

Castro Camilo, D., Lombardo, L., Mai, P.M., Jie, D., Huser, R., 2017. Handling high predictor dimensionality in slope-unit-based landslide susceptibility models through LASSO-penalized generalized linear model. Environ. Model. Softw. 97, 145–156.

Chen, Z., Wang, J., 2007. Landslide hazard mapping using logistic regression model in Mackenzie Valley, Canada. Nat. Hazards 42 (1), 75–89.

Ciurleo, M., Cascini, L., Calvello, M., 2017. A comparison of statistical and deterministic methods for shallow landslide susceptibility zoning in clayey soils. Eng. Geol. 223, 71–81.

Conoscenti, C., Rotigliano, E., Cama, M., Caraballo-Arias, N.A., Lombardo, L., Agnesi, V., 2016. Exploring the eﬀect of absence selection on landslide susceptibility models: a case study in Sicily, Italy. Geomorphology 261, 222–235.

Copas, J.B., Long, T., 1991. Estimating the residual variance in orthogonal regression with variable selection. Journal of the Royal Statistical Society. Series D (The Statistician) 40 (1), 51–59.

Crozier, M., 2005. Multiple-occurrence regional landslide events in New Zealand: hazard management issues. Landslides 2 (4), 247–256.

Derksen, S., Keselman, H.J., 1992. Backward, forward and stepwise automated subset selection algorithms: frequency of obtaining authentic and noise variables. Br. J. Math. Stat. Psychol. 45 (2), 265–282.

Eeckhaut, M.V.D., Marre, A., Poesen, J., 2010. Comparison of two landslide susceptibility assessments in the Champagne–Ardenne region (France). Geomorphology 115 (1), 141–155.

Erener, A., Mutlu, A., Düzgün, H.S., 2016. A comparative study for landslide susceptibility mapping using GIS-based multi-criteria decision analysis (MCDA), logistic regression (LR) and association rule mining (ARM). Eng. Geol. 203, 45–55.

Feng, H., Yu, J., Zheng, J., Tang, X., Peng, C., 2016. Evaluation of diﬀerent models in rainfall-triggered landslide susceptibility mapping: a case study in Chunan, Southeast China. Environmental Earth Sciences 75 (21), 1399.

Frattini, P., Crosta, G., Carrara, A., 2010. Techniques for evaluating the performance of landslide susceptibility models. Eng. Geol. 111 (1), 62–72.

Friedman, J., Hastie, T., Tibshirani, R., 2009. glmnet: LASSO and elastic-net regularized generalized linear models. In: R package version. 1(4).

Goo, Y.-J.J., Chi, D.-J., Shen, Z.-D., 2016. Improving the prediction of going concern of Taiwanese listed companies using a hybrid of LASSO with data mining techniques. Springerplus 5 (1), 539.

Guns, M., Vanacker, V., 2012. Logistic regression applied to natural hazards: rare event logistic regression with replications. Nat. Hazards Earth Syst. Sci. 12 (6), 1937–1947.

Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., Galli, M., 2006. Estimating the quality of landslide susceptibility models. Geomorphology 81 (1), 166–184.

Hansen, A., 1984. Landslide hazard analysis. In: Slope Instability. Wiley, New York, pp. 523–602.

Harrell Jr., F.E., 2015. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. Springer.

Heckmann, T., Gegg, K., Gegg, A., Becht, M., 2014. Sample size matters: investigating the eﬀect of sample size on a logistic regression susceptibility model for debris ﬂows. Nat. Hazards Earth Syst. Sci. 14 (2), 259.

Heerdegen, R.G., Beran, M.A., 1982. Quantifying source areas through land surface curvature and shape. J. Hydrol. 57 (3–4), 359–373.

Hussin, H.Y., Zumpano, V., Reichenbach, P., Sterlacchini, S., Micu, M., van Westen, C., Bălteanu, D., 2016. Diﬀerent landslide sampling strategies in a grid-based bi-variate statistical susceptibility model. Geomorphology 253, 508–523.

Ko, F.W., Lo, F.L., 2018. From landslide susceptibility to landslide frequency: a territory-wide study in Hong Kong. Eng. Geol. 242, 12–22.

Lee, C.-T., 2014. Statistical seismic landslide hazard analysis: an example from Taiwan. Eng. Geol. 182, 201–212.

Lentini, F., Catalano, S., Carbone, S., 2000. Carta geologica della provincia di Messina, Provincia Regionale di Messina. In: Assessorato Territorio–Servizio geologico, SELCA, Firenze.

Leoni, G., Barchiesi, F., Catallo, F., Dramis, F., Fubelli, G., Lucifora, S., Mattei, M., Pezzo, G., Puglisi, C., 2009. GIS methodology to assess landslide susceptibility: application to a river catchment of Central Italy. Journal of maps 5 (1), 87–93.

Lin, G.-F., Chang, M.-J., Huang, Y.-C., Ho, J.-Y., 2017. Assessment of susceptibility to rainfall-induced landslides using improved self-organizing linear output map, support

vector machine, and logistic regression. Engineering Geology224 (Supplement C), 62–74.

Lombardo, L., Cama, M., Maerker, M., Rotigliano, E., 2014. A test of transferability for landslides susceptibility models under extreme climatic events: application to the Messina 2009 disaster. Nat. Hazards 74 (3), 1951–1989.

Lombardo, L., Cama, M., Conoscenti, C., Märker, M., Rotigliano, E., 2015. Binary logistic regression versus stochastic gradient boosted decision trees in assessing landslide susceptibility for multiple-occurring landslide events: application to the 2009 storm event in Messina (Sicily, southern Italy). Nat. Hazards 79 (3), 1621–1648.

Lombardo, L., Bachofer, F., Cama, M., Märker, M., Rotigliano, E., 2016a. Exploiting maximum entropy method and ASTER data for assessing debrisﬂow and debris slide susceptibility for the Giampilieri catchment (North-Eastern Sicily, Italy). Earth Surf. Process. Landf. 41 (12), 1776–1789.

Lombardo, L., Fubelli, G., Amato, G., Bonasera, M., 2016b. Presence-only approach to assess landslide triggering-thickness susceptibility: a test for the Mili catchment (North-Eastern Sicily, Italy). Nat. Hazards 84 (1), 565–588.

Lombardo, L., Opitz, T., Huser, R., 2018a. Point process-based modeling of multiple debrisﬂow landslides using INLA: an application to the 2009 Messina disaster. Stoch. Env. Res. Risk A. 32 (7), 2179–2198.

Lombardo, L., Saia, S., Schillaci, C., Mai, P.M., Huser, R., 2018b. Modeling soil organic carbon with quantile regression: dissecting predictors' eﬀects on carbon stocks. Geoderma 318, 148–159.

Marchesini, I., Ardizzone, F., Alvioli, M., Rossi, M., Guzzetti, F., 2014. Non-susceptible landslide areas in Italy and in the Mediterranean region. Nat. Hazards Earth Syst. Sci. 14 (8), 2215–2231.

Miller, A., 2002. Subset selection in regression. In: Chapman & Hall/CRC Monographs on Statistics & Applied Probability.

Moore, I.D., Grayson, R., Ladson, A., 1991. Digital terrain modelling: a review of hy-drological, geomorphological, and biological applications. Hydrol. Process. 5 (1), 3–30.

Parker, R.N., Rosser, N.J., Hales, T.C., 2017. Spatial prediction of earthquake-induced landslide probability. Natural Hazards and Earth System Sciences Discussions 2017, 1–29.

Petschko, H., Brenning, A., Bell, R., Goetz, J., Glade, T., 2014. Assessing the quality of landslide susceptibility maps–case study Lower Austria. Nat. Hazards Earth Syst. Sci. 14 (1), 95–118.

Pourghasemi, H.R., Rossi, M., 2016. Landslide susceptibility modeling in a landslide prone area in Mazandarn Province, north of Iran: a comparison between GLM, GAM, MARS, and M-AHP methods. Theoretical and Applied Climatology 1–25.

Raja, N.B., Çiçek, I., Türkoğlu, N., Aydin, O., Kawasaki, A., 2017. Landslide susceptibility mapping of the Sera River Basin using logistic regression model. Nat. Hazards 85 (3), 1323–1346.

Rossi, M., Reichenbach, P., 2016. LAND-SE: a software for statistically based landslide susceptibility zonation, version 1.0. Geoscientiﬁc Model Development 9 (10), 3533.

Rossi, M., Guzzetti, F., Reichenbach, P., Mondini, A.C., Peruccacci, S., 2010. Optimal landslide susceptibility zonation based on multiple forecasts. Geomorphology 114 (3), 129–142.

Rotigliano, E., Agnesi, V., Cappadonia, C., Conoscenti, C., 2011. The role of the diagnostic areas in the assessment of landslide susceptibility models: a test in the Sicilian chain. Nat. Hazards 58 (3), 981–999.

Rouse Jr, J., Haas, R., Schell, J. and Deering, D. (1974) Monitoring vegetation systems in the Great Plains with ERTS.

Schillaci, C., Lombardo, L., Saia, S., Fantappiè, M., Märker, M., Acutis, M., 2017. Modelling the topsoil carbon stock of agricultural lands with the Stochastic Gradient Treeboost in a semi-arid Mediterranean region. Geoderma 286, 35–45.

Shou, K.-J., Lin, J.-F., 2016. Multi-scale landslide susceptibility analysis along a mountain highway in Central Taiwan. Eng. Geol. 212, 120–135.

Steger, S., Brenning, A., Bell, R., Glade, T., 2016. The propagation of inventory-based positional errors into statistical landslide susceptibility models. Nat. Hazards Earth Syst. Sci. 16 (12), 2729–2745.

Süzen, M.L., Doyuran, V., 2004. A comparison of the GIS based landslide susceptibility assessment methods: multivariate versus bivariate. Environ. Geol. 45 (5), 665–679.

Tibshirani, R., 1996. Reegression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society. Series B (Methodological) 267–288.

Tziritis, E., Lombardo, L., 2016. Estimation of intrinsic aquifer vulnerability with index-overlay and statistical methods: the case of eastern Kopaida, central Greece. Applied Water Science 1–15.

Van Den Eeckhaut, M., Reichenbach, P., Guzzetti, F., Rossi, M., Poesen, J., 2009. Combined landslide inventory and susceptibility assessment based on diﬀerent mapping units: an example from the Flemish Ardennes, Belgium. Nat. Hazards Earth Syst. Sci. 9 (2), 507–521.

Wang, L.-J., Sawada, K., Moriguchi, S., 2013. Landslide susceptibility analysis with lo-gistic regression model based on FCM sampling strategy. Computers & Geosciences 57 (Supplement C), 81–92.

Zevenbergen, L.W., Thorne, C.R., 1987. Quantitative analysis of land surface topography. Earth Surf. Process. Landf. 12 (1), 47–56.