## 2 Materials and Methods
### Statistical Model, Validation Tools and Model-Building Strategies  
The Multivariate Adaptive Regression Splines (MARS; [43]) method was employed to detect the relations between the outcome and predictors. MARS is based on non-parametric regression, which exploits the splitting of each independent variable into hinge functions to boost the maximum-likelihood-based adaptation skill of the logistic regression method, according to  
\[y=f(x)=\alpha+\sum_{i=1}^{N}\beta_{i}h_{i}(x) \tag{1}\]  
where \(y\) is the dependent variable (the outcome) predicted by the function _f(x)_, \(\alpha\) is the model intercept, and \(\beta_{i}\) is the coefficient of the \(h_{i}\) basis functions, given the \(N\) number of basis functions.  
The MARS method is largely employed in several fields of science (e.g., [44; 45]), but in the last ten years, MARS has been successfully employed in landslide susceptibility evaluation (e.g., [29; 30; 36; 40; 42; 46; 47; 48; 49; 50]). In this research, MARS analysis was performed by exploiting the \"earth\" R-package [51]. Based on previous MARS applications [26; 29; 30; 37; 42; 50], in order to avoid overfitted models and obtain results with geomorphological adequacy in terms of evaluation, we accepted the default setting of the \"earth\" R-package (see [51] for details), limiting the number of terms in the range of 20-200 and setting it with no interaction.  
The MARS method is based on a presence-absence approach. For this reason, the random extraction of the same number of negative and positive cases was carried out. To test the independence of the results (resolution, precision and robustness) from the choice of the selected cases [29; 30; 42], the extraction of the balanced datasets was repeated fifty times. Moreover, a balanced random partition of each dataset into calibration (75%) and validation (25%) subsets was replicated one hundred times. The prediction skill and the reliability of the models was estimated on the base of five thousand validation tests.  
According to [PERSON] and [PERSON] [52], the AUC value (Area under Curve) in the ROC (Receiver Operating Characteristics [53; 54; 55]) was employed to evaluate the prediction skill of the model. At the same time, the Youden index optimised score cut-off [56] was obtained from the ROC plots. In this way, True Positive (TP), True Negative (TN), False Positive (FP) and False Negative (FN) cases were defined, and the related cut-off-dependent indices were calculated. For the balanced analysed datasets, the sensitivity (TPrate) is equal to the ratio of TP/positive cases and 1-FN; the specificity (TNrate) is equal to the ratio of TN/negative cases and 1-FP; the accuracy was equal to the ratio of True/False rate.  
In this research, using five different cell sizes (5, 8, 10, 16 and 32 m) and three types of landslides, a set of 15 models was evaluated. For each model, the prediction skill in evaluating the blinded 25% (prediction), and 100% (success) of the archive was analysed by exploiting all of the validation indices/instruments (ROC plot, AUC value and confusion matrix). On the other hand, an in-depth comparison of the employed geo-environmental variables (or classes) for the construction of the 15 models was carried out. In fact, using the \"varImp\" function of the Caret R-package [57], the percentage of use of each variable by the replicates was calculated.