Comparison of Statistical Analysis Models for Susceptibility Assessment of Earthquake-Triggered Landslides: A Case Study from 2015 Earthquake in Lefkada Island

The main purpose of this study is to comparatively assess the susceptibility of earthquake-triggered landslides in the island of Lefkada (Ionian Islands, Greece) using two different statistical analysis models, a bivariate model represented by frequency ratio (FR), and a multivariate model represented by logistic regression (LR). For the implementation of the models, the relationship between geo-environmental factors contributing to landslides and documented events related to the 17th November 2015 earthquake was investigated by geographic information systems (GIS)-based analysis. A landslide inventory with events attributed to the specific earthquake was prepared using satellite imagery interpretation and field surveys. Eight factors: Elevation, slope angle, slope aspect, distance to main road network, distance to faults, land cover, geology, and peak ground acceleration (PGA), were considered and used as thematic data layers. The prediction capability of the models and the accuracy of the resulting susceptibility maps were tested by a standard validation method, the receiver operator characteristic (ROC) analysis. Based on the validation results, the output map with the highest reliability could potentially constitute an ideal basis for use within regional spatial planning as well as for the organization of emergency actions by local authorities.


Introduction
Globally, landslides are related to a great number of human losses and economic damages ( Table 1). One of the major mechanisms by which a landslide can occur is triggering due to earthquakes [1]. Over the last two decades, strong earthquakes have triggered the occurrence of landslides in many islands of the Ionian Sea (2003 in Lefkada, 2014 in Kephalonia, 2015 in Lefkada, 2018 in Zakynthos). This part of Greece can be characterized as a high-seismicity region, particularly prone and vulnerable to earthquake-triggered landslides.
In order to determine the most fragile landslide areas under the influence of a given earthquake in the future, it is important to indicate these areas through landslide susceptibility assessment and mapping [2]. A landslide susceptibility map gives an indication of "where" future landslides are likely to occur over a region on the basis of local geo-environmental conditions. Susceptibility mapping of earthquake-triggered landslides can be particularly useful for decision makers and planners. The knowledge about the spatial probability of their occurrence is very important for hazard management and mitigation, and safe planning [3].

Study Area
Lefkada is located in the Ionian Sea, between Corfu and Kefalonia, and a narrow strip of sea separates it from mainland Greece (Figure 1). With an extent of about 300 km 2 , it is the fourth largest Ionian Island. Lefkada can be characterized as a mostly mountainous region (about 70% of the total area) with an average altitude of 500 m. Its climate includes mild and rainy winter, and hot summer. Because of these climatic conditions and its geomorphology, the island offers a variety of landscapes, a key element of which is the rich vegetation, narrow plateaus, white beaches, springs, and gorges.

Data
Based on the above-described study area, the appropriate data were collected. Specifically, in order to carry out the landslide susceptibility assessment, a spatial database was designed and developed in GIS environment (by ArcGIS ver. 10.2 software package, ESRI, California, CA, USA). This database included two different types of datasets: (a) The landslide inventory dataset, and (b) the datasets with the geo-environmental factors.

Landslide Inventory
The landslide inventory contains critical information in order to apply quantitative modeling. Following the introduced landslide classification by [29], in the present inventory the term landslide is used to describe earth and rock slides, debris flows, and complex movements (Figure 1). Its creation was based on the use of cloud-free satellite imagery from Google Earth, and field surveys. A Geologically, most of its terrain belongs to the western margin of the Ionian zone covered by thick limestone formations, while a small part of its south-western territory belongs to the northern extent of the pre-Apulian (Paxos) zone composed of limestones extending to a total thickness of 200 m. Even younger deposits have been detected on the north-east of the island, while the most recent deposits comprising alluvial sediments of loose sand and soft marine clays are found in the homonym capital town [28].
In regard to its tectonic aspect, a large number of faults are situated in Lefkada as a result of its geodynamic complexity. However, the main active tectonic structure accommodating the relative motion of the African and Eurasian lithospheric plates, and thus highly influencing the seismicity of the region, is the Kefalonia transform fault (KTF) (Figure 1). The fault consists of two segments, the Kefalonia and Lefkada segments. Its length and width for the Kefalonia segment is estimated as 85 and 20 km, respectively, whereas for the Lefkada segment as 40 and 15 km, respectively. Due to the activity of KTF, Lefkada has a rich history of strong earthquakes with magnitudes up to 6.5 (1914, 1948, 1973, and 2003). An earthquake of similar magnitude occurred in the island on 17th November, 2015.

Data
Based on the above-described study area, the appropriate data were collected. Specifically, in order to carry out the landslide susceptibility assessment, a spatial database was designed and developed in GIS environment (by ArcGIS ver. 10.2 software package, ESRI, California, CA, USA). This database included two different types of datasets: (a) The landslide inventory dataset, and (b) the datasets with the geo-environmental factors.

Landslide Inventory
The landslide inventory contains critical information in order to apply quantitative modeling. Following the introduced landslide classification by [29], in the present inventory the term landslide is used to describe earth and rock slides, debris flows, and complex movements ( Figure 1). Its creation was based on the use of cloud-free satellite imagery from Google Earth, and field surveys. A landslide database maintained by the Laboratory of Engineering Geology in the Department of Geology at the University of Patras [30] was previously exploited for the detection of landslide events triggered by the 17th November 2015 earthquake in the island of Lefkada. The final landslide inventory contained 43 landslide events. Given that in landslide inventories it is not always possible to differentiate their depletion and accumulation zones [15], these zones were mapped together in an entire area forming a single polygon feature for each landslide event ( Figure 1).

Geo-Environmental Factors
A landslide constitutes the result of interaction between several geo-environmental factors. These factors can be subdivided into two categories: (a) The preconditioning factors that are expected to have an effect on the landslide occurrence, and (b) the triggering factors that trigger it. Due to the "nature" of the examined landslides as earthquake-triggered events, both preconditioning and triggering factors were used in this study. The selection of these factors mainly depends on the scale of the analysis and the availability of data [31]. Seven preconditioning factors representing the conditions in the study area before the occurrence of the earthquake, such as land cover, lithology, elevation, slope angle, slope aspect, distance to main road network, and distance to faults, and one (during the) earthquake-related triggering factor like peak ground acceleration (PGA) were taken into account.
Changes in land cover as a result of human activities like forest logging, deforestation, cultivation on steep slopes, and road construction can have a significant impact on landslide occurrence [32]. Land cover types were defined from the European CORINE (Coordinate of Information on the Environment) program [33] by following the level-3 classification scheme of the CORINE 2018 data. Given that the various lithological formations have different slope stability performance, lithology is also characterized as a significant factor for landslide susceptibility assessment [34]. The lithological units of the study area were derived from the geological map sheets (Leykas and Agios Petros) at 1:50,000 scale, provided by the Institute of Geology and Mineral Exploration (IGME). It is worth mentioning that by grouping the initial categories based on their common characteristics, the final land cover and lithology datasets were created.
The elevation factor is useful to classify the local relief and locate points of maximum and minimum heights within terrains [31]. Furthermore, an increase in slope angle is correlated with an increased likelihood of failure [35]. The slope aspect expressing the azimuthal orientation of a slope constitutes a decisive proxy for slope-influencing natural processes like incoming precipitation, differential weathering, as well as wind and solar radiation, which may modulate slope stability. The elevation, slope angle, and slope aspect datasets were produced from the ASTER DEM (digital elevation model) product (with spatial resolution of 30 m) using GIS-based spatial analysis tools.
Road openings at the slope bases have negative impacts on slope stability. The dataset of distance to main road network was based on the digitization of the main roads from satellite images. Moreover, since the seismic ground acceleration and intensity diminish with the distance to tectonic elements, it was necessary for the relative factors (distance to faults and PGA) to be included in the present study. The dataset of distance to faults was derived from the digitization of the faults from the 1:50,000 geological map sheets. A PGA dataset was obtained from the earthquake strong-motion record of the Hellenic Strong Motion Network (HSMN) maintained by the Institute of Geodynamics-National Observatory of Athens (IG-NOA). At a national scale, HSMN consists of permanent stations in continuous and triggering mode across the entirety of Greece, an instrument-based observation array, and a network management system [36].
All the factors were organized in the relative GIS data layers and converted into raster grid format with pixel size 30 m. An overview of these factors is presented in Figure 2. depends on the scale of the analysis and the availability of data [31]. Seven preconditioning factors representing the conditions in the study area before the occurrence of the earthquake, such as land cover, lithology, elevation, slope angle, slope aspect, distance to main road network, and distance to faults, and one (during the) earthquake-related triggering factor like peak ground acceleration (PGA) were taken into account. Changes in land cover as a result of human activities like forest logging, deforestation, cultivation on steep slopes, and road construction can have a significant impact on landslide occurrence [32]. Land cover types were defined from the European CORINE (Coordinate of Information on the Environment) program [33] by following the level-3 classification scheme of the CORINE 2018 data. Given that the various lithological formations have different slope stability performance, lithology is also characterized as a significant factor for landslide susceptibility assessment [34]. The lithological units of the study area were derived from the geological map sheets (Leykas and Agios Petros) at 1:50,000 scale, provided by the Institute of Geology and Mineral Exploration (IGME). It is worth mentioning that by grouping the initial categories based on their common characteristics, the final land cover and lithology datasets were created. The elevation factor is useful to classify the local relief and locate points of maximum and minimum heights within terrains [31]. Furthermore, an increase in slope angle is correlated with an increased likelihood of failure [35]. The slope aspect expressing the azimuthal orientation of a slope constitutes a decisive proxy for slope-influencing natural processes like incoming precipitation, differential weathering, as well as wind and solar radiation, which may modulate slope stability. The elevation, slope angle, and slope aspect datasets were produced from the ASTER DEM (digital elevation model) product (with spatial resolution of 30 m) using GIS-based spatial analysis tools.
Road openings at the slope bases have negative impacts on slope stability. The dataset of distance to main road network was based on the digitization of the main roads from satellite images. Moreover, since the seismic ground acceleration and intensity diminish with the distance to tectonic elements, it was necessary for the relative factors (distance to faults and PGA) to be included in the present study. The dataset of distance to faults was derived from the digitization of the faults from the 1:50,000 geological map sheets. A PGA dataset was obtained from the earthquake strong-motion record of the Hellenic Strong Motion Network (HSMN) maintained by the Institute of Geodynamics-National Observatory of Athens (IG-NOA). At a national scale, HSMN consists of

Frequency Ratio (FR) Model
Frequency ratio is a simple and understandable bivariate statistical model based on the spatial associations between distribution of landslide events and each of the landslide-influencing factors, to expose the level of their correlation. It is calculated for each category of the factors by dividing the landslide occurrence ratio by the area ratio: where N pix (S j ) is the number of landslide pixels in factor category j, and N pix (N j ) is the number of pixels in the same factor category. A FR value of 1 (average value) means that the density of landslides in the category is proportional to the size of the category. If the value is greater than 1, then there is a high correlation, whereas a value of less than 1 means a lower correlation [37].
The overall landslide susceptibility (LS) value for each pixel can be obtained by summing the FR values of different landslide-influencing factors: where FR i,j is the frequency ratio value for the category j of the factor i, and n is the total number of the factors.

Logistic Regression (LR) Model
A logistic regression model forms a multivariate relation between a dependent variable and several independent variables [38]. The dependent variable is binary (i.e. it can take only the value 1 and 0), while its predictors, the independent variables, can be either continuous, discrete, or any combination of both types. In terms of a landslide susceptibility assessment, the goal of LR is to find the best fitting model to describe the relationship between the absence or presence (value of 0 or 1) of landslides (dependent variable) and a set of landslide-influencing factors (independent variables). The model can be expressed in its simplest form as: where P is the probability of landslide occurrence, which ranges from 0 to 1, and z is a linear sum of a constant and the product of the independent variables with their respective coefficients. The value of z varies from -∞ to +∞ and is calculated from the equation: where n is the number of independent variables, x i (i = 1, 2, . . . , n) are the independent variables, b 0 is the constant of the model, and b i (i = 1, 2, . . . , n) are the coefficients. Another way to write the LR model is the logit transformation using an equation in the form of: where p(Y = 1)/1 − p(Y = 1) is the so-called odds or likelihood ratio. The result represents the probability that a landslide event will occur divided by the probability that it fails to do so. If a coefficient is positive, its transformed log value will be greater than 1, meaning that the event is more likely to occur. If a coefficient is negative, the relative value will be less than 1 and the odds of the event occurring decreases. A coefficient of 0 has a transformed log value of 1, and it does not change the odds one way or the other [31]. Generally, the model estimates the coefficients and statistics based on the values of independent variables and the status of the dependent variable in a sampling dataset, using a maximum likelihood method [39]. Using the estimates derived from the implementation of the model on the selected sample, the probability of a landslide can be calculated on a pixel by pixel basis (Equation (3)).

Data Processing
The data processing included two different processes, such as the landslide data sampling and preparation of factor data. For the landslide data sampling, the landslide dataset derived from the inventory should be separated into two parts: (a) A training dataset, for the implementation of the models, and (b) a validation dataset, for the evaluation of their outputs. The total of 43 landslide polygons were converted into raster grid format with pixel size 30 m, resulting to 542 pixels. Among these pixels, 80% (434 pixels) was randomly selected as the training dataset, and the remaining 20% (108 pixels) made up the validation dataset [2,3]. Furthermore, an equal number of pixels from the landslide-not-occurrence area were randomly selected for both the training and validation datasets. The FR model handled only a landslide dataset, whereas the LR model handled both the landslide and non-landslide datasets. Thus, the training dataset for FR model contained 434 pixels, and for the LR model 868 pixels. The validation dataset contained a total of 216 pixels. In both training and validation datasets, the target value of 1 was assigned to the landslide pixels, while the target value of 0 to the non-landslide pixels.
For the preparation of factor data used in FR model, the GIS-based "Natural Breaks" categorization was implemented for most of the factors with a continuous numerical scale (elevation, slope angle, distance to main road network, and distance to faults), except for PGA factor whose categorization was executed in a manually according to its presented values. In "Natural Breaks" categorization, category breaks identify best group similar values and maximize the differences between categories according to the deviations about the median [40]. For the preparation of factor data used in the LR model, the factors with continuous numerical scale, such as elevation, slope angle, distance to main road network, distance to faults, and PGA were handled in their original format, in order not to alter the state and information presented in their data layers. On the contrary, the factors with discrete categorical scale, such as land cover, lithology, and slope aspect were re-scaled in the range of 0.1 to 0.9 by coding and ranking their categories based on the relative landslide densities.

Implementation of Models
For the FR model, by crossing the landslide training dataset (434 pixels) with each factor layer, the relative landslide density in each factor category was calculated. The FR value for each of these categories was then estimated with Equation (1) ( Table 2). Finally, all factor layers were overlaid by summation (Equation (2)) and a resultant the landslide susceptibility map was obtained. This map was categorized into five categories ("Very Low," "Low," "Moderate," "High," and "Very High" susceptibility) through the "Natural Breaks" method ( Figure 3). Table 2. Values derived from the frequency ratio (FR) model for all factor categories, and multicollinearity checking indexes (tolerance (TOL) and variance inflation factor (VIF)) and coefficients derived from logistic regression (LR) model for all factors. and 1% of the ''High'' and ''Very High'', categories of the FR model were characterized with ''Very Low'' susceptibility in the map from the LR model. Concerning the coverage similarities, it is indicated that 2% of the study area was characterized by similar susceptibility for each of two, "High" and "Very High", categories. The corresponding percentages for the other three, "Moderate", "Low", and "Very Low", categories were 3%, 4% and 21%, respectively. The distribution of these similarities are presented in Figure 5.     In case of the LR model, the training dataset (868 landslide and non-landslide pixels) were matched with each factor layer in order to create a database. This database with the eight factors as independent variables, and the presence and absence of landslide (binary value of 0 and 1) as the dependent variable was imported into the SPSS (ver. 23, Company, City, Country IBM, New York, N.Y., USA) software package. The tolerance (TOL) and variance inflation factor (VIF) indexes were then estimated to check the multicollinearity of independent variables. TOL estimation is based on the determination of a regression of each explanatory variable on all the other explanatory variables. VIF is calculated by 1/TOL [41]. These indexes were found to be greater than 0.2 for TOL and less than 10 for VIF revealing that there is no multicollinearity between any of the factors ( Table 2). Following the multicollinearity checking, the binary LR algorithm was carried out to calculate the correlation of landslides to each factor. This correlation is expressed by the coefficients (Table 2). After the assignment of coefficients to all the factors, a GIS-based weighted overlay was applied using Equation (4). By inserting the output of overlay into Equation (3), a resultant landslide susceptibility map was produced and finally categorized into five categories (Figure 3).

Results
The results of the FR and LR models are summarized in Table 2. Based on the FR model, the factor categories that presented the highest correlation with the occurrence of landslides triggered by the 17th November 2015 earthquake were "0.12 to 0.16 g" (FR = 3.47) for PGA, "west" facing (FR = 3.35) for slope aspect, "open spaces with little/no vegetation" (FR = 3.24) for land cover, "alluvium deposits" (FR = 2.89) for lithology, "0 to 134 m" (FR = 2.78) for elevation, ''35 to 65 degrees" (FR = 2.68) for slope angle, and "0 to 324 m" (FR = 2.24) for distance to faults. There is no mention for the factor of distance to main road network because its highest value (FR = 1.66) was very close to 1 (average value). According to the LR model, the factors of PGA, slope aspect, land cover, lithology, and slope angle had positive coefficients, whereas the factors of elevation, distance to main road network, and distance to faults had negative coefficients. Among the factors with positive coefficients, the highest value was presented by PGA (coefficient of 60.946) followed by slope aspect (coefficient of 4.709).
The landslide susceptibility maps produced by the two models are presented in Figure 3. In the map from the FR model, the "High" and "Very High" susceptibility categories were mainly located in the western coastal zone of Lefkada Island, including large pockets of "High" susceptibility in its central and southern parts. These two categories cover 13% and 3%, respectively, of the study area ( Figure 4). In the map from the LR model, the "Very High" susceptibility category was detected in the western coastal zone of the island, with some limited pockets of "High" and "Very High" susceptibility in its southern part. The coverage percentages for the two susceptibility categories were found to be 4% and 5%, respectively ( Figure 4). Moreover, the cross-comparison of the two output maps in terms of coverage differences between their categories (Table 3) indicated that the highest difference percentages (29% and 19%, respectively) were shown between the "Low" and "Moderate" categories of the FR model, and the "Very Low" category of LR model. It is also worth mentioning that no coverage relation was observed between the "High" and "Very High" categories of the LR model, and the "Very Low" and "Low" categories of the FR model. On the contrary, 5% and 1% of the "High" and "Very High", categories of the FR model were characterized with "Very Low" susceptibility in the map from the LR model. Concerning the coverage similarities, it is indicated that 2% of the study area was characterized by similar susceptibility for each of two, "High" and "Very High", categories. The corresponding percentages for the other three, "Moderate", "Low", and "Very Low", categories were 3%, 4% and 21%, respectively. The distribution of these similarities are presented in Figure 5.   Table 3. Coverage cross-comparison for the landslide susceptibility categories between both the FR and LR models.

Validation of Results
Validation is an essential process to obtain knowledge about the accuracy and prediction ability of the models [42]. The overlay of produced susceptibility maps with a landslide dataset is considered as a standard validation method. In the present study, the results of this overlay with the landslide validation dataset ( Figure 6) showed that, for the FR model, 69% and 28% of the landslide pixels were found within the ''Very High'' and ''High'', susceptibility categories and none of these pixels were found within the ''Very Low'' and ''Low'' susceptibility categories. Regarding the LR model, the percentages were 87% and 8% for the ''Very High'' and ''High'', susceptibility categories, and only 1% for the ''Very Low'' and ''Low'' susceptibility categories.
A more advanced validation method, which has been widely applied, is the receiver operating characteristics (ROC) analysis [8,15,43,44]. Based on the overlay of landslide susceptibility maps with an "independent" validation dataset (with landslide and non-landslide data), a ROC graph was created. In this graph, the sensitivity of a model, which is determined as the percentage of the correctly predicted landslide pixels by the model, was plotted against specificity, which is the percentage of predicted landslide pixels over the total study area [45]. The value of area under the ROC curve (AUC) indicates the prediction ability of the model. With a range from 0.5 to 1.0, the higher this value is, the better is the prediction ability of the model. More details about ROC analysis can be found in [46]. The results of ROC analysis are presented in Figure 7 and Table 4. The AUC value of the LR model was equal to 0.984, followed by the FR model with a value of 0.976. The accuracy percentages were equal to 93.1% and 78.2%, respectively.

Validation of Results
Validation is an essential process to obtain knowledge about the accuracy and prediction ability of the models [42]. The overlay of produced susceptibility maps with a landslide dataset is considered as a standard validation method. In the present study, the results of this overlay with the landslide validation dataset ( Figure 6) showed that, for the FR model, 69% and 28% of the landslide pixels were found within the "Very High" and "High", susceptibility categories and none of these pixels were found within the "Very Low" and "Low" susceptibility categories. Regarding the LR model, the percentages were 87% and 8% for the "Very High" and "High", susceptibility categories, and only 1% for the "Very Low" and "Low" susceptibility categories.
A more advanced validation method, which has been widely applied, is the receiver operating characteristics (ROC) analysis [8,15,43,44]. Based on the overlay of landslide susceptibility maps with an "independent" validation dataset (with landslide and non-landslide data), a ROC graph was created. In this graph, the sensitivity of a model, which is determined as the percentage of the correctly predicted landslide pixels by the model, was plotted against specificity, which is the percentage of predicted landslide pixels over the total study area [45]. The value of area under the ROC curve (AUC) indicates the prediction ability of the model. With a range from 0.5 to 1.0, the higher this value is, the better is the prediction ability of the model. More details about ROC analysis can be found in [46]. The results of ROC analysis are presented in Figure 7 and Table 4. The AUC value of the LR model was equal to 0.984, followed by the FR model with a value of 0.976. The accuracy percentages were equal to 93.1% and 78.2%, respectively. Geosciences 2019, 9, x FOR PEER REVIEW 13 of 17

Discussion
For the seismically active region of Ionian Islands, and specifically the island of Lefkada, the occurrence of strong earthquakes acting as a triggering mechanism for landslides, is a frequent phenomenon. Therefore, the acquisition of knowledge about the spatial definition of potential earthquake-triggered landslides is considered as an essential need for the region. In terms of acquisition of this knowledge, the present study focused on the creation of a reliable landslide susceptibility map that will spatially define the landslide occurrence for the island of Lefkada. This was designed in case an earthquake with a magnitude similar to the one witnessed on the 17th November 2015 occurs in the future.
In order to obtain the desired map, two different statistical models, such as bivariate FR and multivariate LR were applied based on the integration of eight geo-environmental factors and a landslide inventory including landslides triggered by the specific earthquake. As a bivariate model, FR uses the observed densities of landslide pixels within factor datasets to estimate susceptibility in a "sub-factor", category level. Concerning the category results of the FR model in this study, the areas of island with low elevation, very high slope angle and western facing, being very close to faults, covered by little or no vegetation and alluvium deposits, and presenting high PGA during the earthquake, were highly correlated with the occurrence of the resultant landslides.
Multivariate LR analyzes variables that are non-symmetrical and show skewed distributions to create a mathematical model that predicts the probability of occurrence in a sampling unit. The coefficients assigned to influencing factors from the LR model are useful to assess their importance on the presence or absence of landslides. In the present study, the factors of PGA, slope aspect, land cover, lithology, and slope angle were found to have positive coefficients indicating that these factors are positively related with the occurrence of earthquake-triggered landslides. However, among them, the factor with the significantly strongest effect on their occurrence was the PGA. It could be said that this finding confirms the "nature" of landslides under investigation as earthquake-triggered events.
The output maps of the two models visualized the spatial distribution of the estimated landslide susceptibility in the study area. Both maps presented the western coastal zone as the part of the island with the highest probability to produce landslides under the influence of a potential earthquake with a magnitude similar to earthquake examined in this paper. According to the study [23], the same part of the island was proved to be riskiest area for earthquake-triggered landslides after the 2003 earthquake.
In regard to reliable results, the susceptibility maps produced by the two models were compared by means of validation using a dataset not originally included in the training process. From the method of overlay, it was derived that the "Very High" and "High" susceptibility categories, contained a total of about 95% of landslide data included in the validation dataset. Although both models present a very good prediction ability (based on their AUC values in Figure 7), the better one is found in the the LR model. The same model also shows superiority against the FR model in terms of accuracy, as the last one tended to present a susceptibility overestimation, which is expressed by high sensitivity and simultaneously low specificity values (Table 4). Both models can successfully detect the susceptible areas ("High" and "Very High" categories), but only the LR model can, at the same time, successfully detect the non-susceptible areas ("Very Low" and "Low" categories). Thus, LR can be assumed as more balanced model against FR. In general, among statistical analysis models, multivariate LR has been proved to give better validation results compared to bivariate models for landslide susceptibility assessments [47][48][49]. Regarding specifically the susceptibility assessment of earthquake-triggered landslides, in the study [2], LR model presented very good validation results from two susceptibility maps produced for a region of China using multi-temporal and earthquake triggered landslide datasets. Moreover, in the study [3], the susceptibility of earthquake-triggered landslides was comparatively evaluated for a watershed of China based on the implementation of six different quantitative models. One of these models was the LR, which provided the best validation results when compared to the other models.
It is worth mentioning that in this study important analysis parameters, such as the number and type of geo-environmental factors, their "pre-model" processing, the size of training and validation datasets, and the categorization of final susceptibility maps, were based on objective criteria. Considering the data-driven "nature" of the two models, the examination of alternative choices for these parameters could lead to different results.

Conclusions and Outlook
Due to its confirmed higher reliability, the landslide susceptibility map produced by the LR model could constitute an important basis for the assessment and management of landslide hazard over the island of Lefkada. This map provides useful information for decision makers and planners to choose suitable locations to implement reconstructions and developments. Furthermore, it could be used for the estimation of the amount of damages and losses that landslide hazards may potentially cause, and the consequent planning of mitigation actions that will reduce the risk of structural damage and loss of life on the island.
Future research work may include testing the LR model in other areas affected by earthquake-triggered landslide events, and its comparison with other advanced soft computing models, such as artificial neural networks (ANNs), support vector machines (SVMs), as well as ensemble random forest and neuro-fuzzy models.