Spatial Analysis on the Variances of Landslide Factors Using Geographically Weighted Logistic Regression in Penang Island, Malaysia

Abstract

Landslides are one of the common natural disasters involving mostly movement of soil surfaces associated with gravitational attraction. Their adverse losses and significant damage, which always result in at least 17% of casualties and billions of dollars of financial losses worldwide, have made landslides the third most notorious phenomenon devastating many parts of the world. Malaysia has had multiple landslide occurrences, particularly in highly urbanized areas, such as Penang Island, owing to the declining vegetation cover in hilly terrains. Thus, this study aims to delineate the spatial relationship variances between landslide occurrences and the influencing factors in the area of interest. Ten influencing factors considered, including distance to roads, distance to rivers, distance to faults, slope angle, slope aspect, curvature, rainfall annual average, lithology, soil series, and land use. In this study, we use a software (GWR 4.0) as a medium for the analysis processing, coupled with GIS. A local statistical technique, Geographically Weighted Logistic Regression (GWLR), is primacy in capturing the geographical variation of the model coefficients that considers non-stationary variables and models their relationships, as well as processes regression coefficients over space. Goodness-of-fit criteria were used to evaluate the GWLR model, namely AICc that decrease from 872.202167 to 800.856998. Bayesian Information Criterion (BIC) shows a decrease in value from 925.784185 to 945.196942. Likewise, deviance decreased from 849.931675 to 739.175630, while pdev increased from 0.379457 to 0.460321. These goodness-of-fit criteria values express GWLR as the best model for local measure. The variances in both local parameter estimates and the t-values (negative and positive values) show the level of significance for each landslide factor in influencing landslide occurrences across the study area. The results of the local parameter estimates and the t-values also show that the slope angle and the slope aspect spatially affect landslide occurrences across the study area. Therefore, a proper perspective and a thorough understanding of the certain slope condition must be established for future mitigation actions to support the agenda of SDG 15, which promotes resilience and disaster risk reduction.

1. Introduction

A landslide is a natural disaster involving mostly the movement of soil surfaces, rock falls, or a combination of both movements, and is commonly associated with gravitational attraction [1,2,3,4]. The number of fatalities often depends on the severity of the landslide event, concerning thousands of human lives and the economy. Although landslides have caused a lower death toll compared to other disasters, the destructive effects of the later consequences are rather devastating. Due to their adverse losses and significant damage, which always result in at least 17% of casualties and billions of dollars of financial losses worldwide, landslides are quickly acknowledged as the third most notorious phenomenon devastating many parts of the world [5,6]. For instance, the mountainous areas of Peru experienced the most catastrophic landslides in 1970, when these events killed approximately 18,000–20,000 people [7]. According to NASA, more than 8935 landslides took place around the world, with 1120 cases scattered in Southeast Asia. On top of that, in the year 2010, the highest number of people affected by landslides was recorded, 89% of which were from China (Figure 1).

Being one of the developing countries in Southeast Asia, Malaysia has had multiple landslide occurrences, but particularly in highly urbanized areas such as Kuala Lumpur, Selangor, and Penang Island. This may be due to the landscape continuously shifting, which mostly involves human activities [9,10,11]. According to [12], the level of urbanization in Penang Island is approximately 90.8% and is among the highest, after Kuala Lumpur, Putrajaya, and Selangor. The growth in urbanization resulted in Penang Island being the most populated island in Malaysia. These days, the significant rise in population has indirectly limited the amount of flat land, thus leading to hillside development. The unrestricted development is clearly agonizing, as areas often undergo slope cutting and land clearing for housing, hotels, and apartments [13]. However, different factors can influence landslides and occasionally correspond to one another. Therefore, this study aims to delineate the spatial relationship between landslide occurrences and the influencing factors in the area of interest, Penang Island.

Comprehensive analysis of landslides can be carried out by qualitative and quantitative methods. The qualitative method is a relatively subjective approach that is entirely based on expert knowledge and opinions and uses the heuristic method technique [14,15,16]. On the other hand, quantitative approaches are mostly based on statistical techniques that can be classified as one of two sub-techniques: bivariate or multivariate analysis. Bivariate analysis is the simplest form of statistical analysis, which is commonly used when two variables are observed against each other; one variable is dependent and the other is independent [17,18]. Bivariate analyses such as Weighted Overlay, Weights of Evidence (WoE), and Frequency Ratio (FR) are widely used in landslide studies [19,20,21,22,23]. Meanwhile, multivariate analysis techniques are purposely employed to observe the relationship between the dependent variable and several independent variables. Logistic Regression (LR), Artificial Neural Networks (ANN), and Analytical Hierarchy Process (AHP) are among the most prevalent multivariate analyses in landslide fields [24,25,26,27,28,29].

In this study, Geographically Weighted Logistic Regression (GWLR), coupled with Geography Information System (GIS) software, is used to analyze 10 influencing factors, taking into consideration landslide spatial analysis. GWLR is an extended version of the standard LR, which can be fully utilized on a local scale. Moreover, GWLR has a significant advantage in binary classification, which is fitting, as this study uses the values 0 and 1 to represent the existence of landslide occurrences. In addition, each landslide occurrence is closely associated with the characteristics of the environment, thus the contribution from factors is less likely identical across the whole of Penang Island. GWLR has also been widely used in landslide studies, which resulted in the best accuracy compared to other models. Therefore, using a heterogeneously designed spatial method, such as GWLR, will generate the best-fitting model for each contributed factor in the study area.

2. Materials and Methods

2.1. Study Area

Penang Island was chosen as the study area due to its history of frequent landslide occurrences. Penang is located in Northwest Peninsular Malaysia at 5.4164° N, 100.3327° E and is 1 of the 14 states of Malaysia, see Figure 2. The terrain consists of coastal plains, hills, and mountains. Penang Island covers an area of approximately 293 km², the fourth largest island in Malaysia, and is separated from the mainland by a channel. It has seen rapid industrialization, rapid technological progress, and the greatest rate of urbanization, resulting in increased demand for land [30]. The Institute of Strategic and International Studies and Penang Development Corporation stated that Penang Island is made up of steep topography above 60 m (exceeding 50% of the land area), which makes half of Penang’s land vulnerable to landslides. Since September 1995, Penang Island has recorded 60 landslide incidents in the Penang Hill area following a freak storm that damaged the pathway near the Penang Botanical Gardens [31]. The changes in land use in Penang Island, which developed along the hillside to accommodate the population growth, are indirectly increasing the frequency of landslide occurrences in the area. This, in the end, will cause an expansion of urbanization. Along with Kuala Lumpur, Penang Island has also experienced rapid urbanization due to industrial activities. For instance, according to the Department of Statistics Malaysia, the level of urbanization in Penang is 90.8% and is among the highest, after Kuala Lumpur, Putrajaya, and Selangor. The saturated urbanization in Penang today is associated with the drastic population increase. Additionally, [32] stated that with the increase in the human population, people have started to dig deeper pits to seek essential minerals. Hence, landslides will be intensified by abrupt and large-scale changes in land use, such as developing new land for neighborhoods, infrastructure, and cities in order to support population demands.

2.2. Flowchart

Figure 3 shows the flowchart of GWLR analysis implemented in the GWR 4.0 software.

2.3. Data Acquisition

The landslide inventory for this study was created based on the databases from previous landslide studies, the database managed by The Public Works Department Malaysia, and field surveys. A total of 988 landslide points (494 landslides, 494 non-landslides) were gathered and used to form this inventory. The dependent variable is binary and is classified with the value 1 (landslide presence) or 0 (landslide absence). The landslide absence points were created randomly using ArcGIS. According to [33], the majority of the landslides that occurred in the center of the island within steep terrain are primarily composed of shallow rotational debris slides and debris flows. Moreover, [17] stated that hill land development in the study area was responsible for most of the landslides over the decades.

2.4. Landslide Influencing Factors

This study collected data representing factors affecting a landslide from different sources (

Table 1. Data sources.

Data	Sources
Landslides	Public Works Department and field survey
Road network	Open source
River network	Department of Irrigation and Drainage
Contour map	Department of Survey and Mapping Malaysia (JUPEM)
Land use map	PLANMalaysia
Geology map	Department of Minerals and Geoscience Malaysia
Soil Series map	Department of Agriculture
Rainfall daily measurements	Department of Irrigation and Drainage

). One study [16] stated that there are no specific rules for choosing landslide characteristics. Thus, a total of 10 influencing factors were analyzed and generated for data–layer maps using GIS software, including distance to roads, distance to rivers, distance to faults, slope angle, slope aspect, curvature, rainfall annual average, lithology, soil series, and land use (Figure 4). Slopes, aspects, and curvatures were derived from a 30 × 30 digital elevation model (DEM) using a contour map.

2.5. Multicollinearity Analysis

The premise of establishing a regression model is that each explanatory variable is independent of the others [34]. This analysis requires the consideration of a possible strong linear correlation between the variables, which may later lead to multicollinearity issues. These issues cause instability in the calculation of the variables resulting in irrelevant results. Thus, the multicollinearity between these 10 variables was considered using two indicators, namely tolerance (TOL) and Variance Inflation Factor (VIF). If the value of TOL is <0.10 and VIF is >5, the results would indicate strong multicollinearity among the factors [35,36]. The following equations were used to measure the TOL and VIF:

$$\mathrm{TOL}=1-R_j^2$$

(1)

$$\mathrm{VIF}=\frac{1}{\mathrm{TOL}}$$

(2)

where $R_j^2$ is the coefficient of regression of j on the added conditioning factors.

2.6. Geographically Weighted Logistic Regression (GWLR)

GWLR is a local statistical technique that considers non-stationary variables and models their relationships, which processes regression coefficients over space [23,37,38,39]. In other words, it is an outgrowth version of the standard Logistic Regression (LR), which determines the relationships between the dependent and independent variables down to a local scale instead of a global one. It was initially developed to traverse the relations between riverbank erosion and geomorphological variables [34]. However, GWLR is now gradually being applied to other studies in science [40]. GWLR analysis for this study utilizes the GWR 4.0 software [41]. The software prefers a spatial database, thus it is constructed first within a GIS software, which consists of the inventory of landslides and analyzed factors. The regression equation can be defined as:

$$P_{\left(u_i,v_i\right)}=\frac{\exp \left({\beta }_{0\left(u_i,v_i\right)}+{\beta }_{1\left(u_i,v_i\right) }{X}_1+\cdots +{\beta }_{k\left(u_i,v_i\right) }{X}_k\right)}{1+\exp \left({\beta }_{0\left(u_i,v_i\right)}+{\beta }_{1\left(u_i,v_i\right) }{X}_1+\cdots +{\beta }_{k\left(u_i,v_i\right) }{X}_k\right)}$$

(3)

where P_{(ui, vi)} is the probability of landslide events occurring at location i with coordinates (u_i, v_i); β₁(u_i, v_i) is the local parameter coefficient of the landslide locations, (u_i, v_i); and β₀(u_i, v_i) is the intercept of point (u_i, v_i) [41].

The stages for the creation of a landslide susceptibility map according to the GWLR model defined above are similar to the global model. The main difference is the application of the logistic regression analysis at a local scale. This is performed by fitting a weighted logistic model for each landslide location using a subset of the original dataset. Thus, a kernel that determines the size of the subset and a weighting function to calculate the weights are required in order to define the GWLR model. In this study, a Gaussian adaptive kernel was adopted. The optimal number of nearest neighbors is estimated by the GWR 4.0 software using an optimization criterion (i.e., the minimization of the corrected Akaike Information Criterion (AICc)). The equation is expressed as:

$$W_{ij=}=\exp \left(-\left(\frac{d_{ij}}{b}\right)2\right)$$

(4)

where W_ij is the weight for unit j in the neighborhood of unit i, d_ij is the distance between the center point of unit i and j as the measurement of spatial proximity degree, and b is the bandwidth of the Gaussian kernel function. Many approaches are available for determining the bandwidth, including the cross-validation (CV) method and the Akaike Information Criterion (AIC) method. However, [42] stated that the AICc method is easier to avoid over-fitting compared to the CV method, hence this study used the AICc method to determine the bandwidth.

3. Results

The multicollinearity analysis of each independent variable is important to reduce over-fitting in the models. Results in

Table 2. Multicollinearity diagnostics between explanatory variables.

Explanatory Variables	Collinearity Statistics
	Tolerance	VIF
Distance to road	0.811	1.234
Distance to river	0.715	1.399
Distance to fault	0.720	1.389
Aspect	0.943	1.061
Slope	0.613	1.632
Curvature	0.973	1.028
Land use	0.921	1.086
Lithology	0.603	1.659
Rainfall annual average	0.957	1.045
Soil	0.755	1.324

below show the lowest TOL value is 0.603 and the highest VIF value is 1.659, which are both laid down on lithology. These values achieved critical values (i.e., TOL < 0.10 or VIF > 5), indicating there is no significant multicollinearity among the 10 landslide factors considered in this study.

Four goodness-of-fit criteria provided by the GWR 4.0 software were used to evaluate the GWLR model. The optimal number of nearest neighbors is estimated using an optimization criterion resulting in the minimization of the corrected AICc.

Table 3. Comparison diagnostics of global regression and GWLR models.

Criteria	Models
	Global Regression	GWLR
Deviance	849.931675	739.175630
AICc	872.202167	800.856998
BIC/MDL	925.784185	945.196942
pdev	0.379457	0.460321

summarizes the results of global and GWLR diagnoses. According to [43], if the difference in the AICc index of the GWLR is >4, the model is considered to have improved. The criterion of the AICc for local regression is 872.202167, which has a lower value than the global regression (800.856998), revealing that GWLR is the best model. Furthermore, the percent deviance (pdev) value increased from 0.379457 to 0.460321, showing that GWLR interprets the spatial variability of landslides better than the global parameter estimates.

Meanwhile,

Table 4. Coefficient values for each variable.

Variable	Coefficients
	Mean	STD	Min	Max	Range
Intercept	7.078098	4.332953	1.909403	17.832890	15.923487
DistRoad	−0.003920	0.000640	−0.005496	−0.002598	0.002898
DistRiver	−0.001008	0.000776	−0.003240	0.000315	0.003555
DistFault	−0.000575	0.000249	−0.001467	−0.000168	0.001300
Aspect	0.001856	0.000979	−0.000331	0.003078	0.003409
Slope	0.013108	0.029410	−0.034743	0.054124	0.088867
Curvature	−0.319739	0.167111	−0.730562	−0.040148	0.690413
Rainfall average	−0.000329	0.000652	−0.001670	0.000592	0.002263
Soil series	−0.942703	0.484519	−2.336661	−0.419098	1.917564
Lithology	−0.407150	0.217785	−0.986768	0.115380	1.102149
Land use	−0.070563	0.085388	−0.293905	0.073151	0.367056

expresses the descriptive statistics of local parameter coefficients of the landslide factors. The intercept ranges from 1.91 to 17.83 with a mean of 7.08. For the Euclidean distance factors: distance to roads ranges from −0.005 to −0.002 with a mean of −0.003; distance to rivers ranges from −0.003 to 0.0003 with a mean of −0.001; and distance to faults ranges from −0.001 to 0.001 with a mean of −0.005. Meanwhile, for the DEM-derived factors: slope aspect ranges from −0.0003 to 0.003 with a mean of 0.001; slope angle ranges from −0.03 to 0.05 with a mean of 0.01; and curvature ranges from −0.73 to −0.04 with a mean of −0.32. For one of the hydrological factors, namely rainfall annual average, the results show a range from −0.0017 to 0.00059 with a mean of −0.00033. Soil series and lithology factors show ranges of −2.34 to −0.42 with a mean of −0.94 and −0.99 to 0.12 with a mean of −0.41, respectively. Lastly, the land use factor ranges from −0.29 to 0.073 with a mean of −0.071.

GWLR is highly efficient when there is a spatial correlation between both variables. In most GWLR applications using multiple independent variables, parameter estimates and t-values are used to determine the influential regions, in which, the values detect the significance level of the relationship between the dependent variable and the independent variables [44]. For the selected independent variables in the GWLR model, raster maps were generated from the parameter estimates and t-values using an Inverse Distance Weighted (IDW) interpolation to show how the relationships between the landslide and each factor vary across space.

4. Discussion

The coefficient maps of the influencing factors are the stepping-stone for landslide mitigation in the study area. The results portrayed in Figure 5 and Figure 6 indicate a spatial variation in the relationship between landslides and the influencing factors in the study area. Because a landslide is a complex phenomenon, its occurrence may due to one or more factors that are constantly associated with each other. It is crucial to identify the factors that are most likely to influence one particular event. The results from the GWLR model reflect in more detail the contribution of each factor to the landslides. Because the spatial heterogeneity of the influencing variables is completely taken into account, the spatial relationship between the landslides and the influencing factors is no longer a fixed value.

For instance, Figure 5 expresses the estimated parameter values calculated in the GWR software. Evaluated values for a particular independent variable that are negative indicate a negative impact on the dependent variable, whereas a positive sign shows a significant influence on the occurrences of the dependent variable. Figure 5a represents parameter estimate values for distance to roads, which constantly shows negative values in all parts of the study area, ranging from −0.0054 to −0.0025. This shows that distance to roads does not influence the occurrence of landslides, as road networks are concentrated on the east side of the study area where the slope’s steepness is less than 15°.

Distance to rivers values, as shown in Figure 5b, convey various values ranging from −0.0032 to 0.00031, with negative values scattered across Penang Island, except for in the southwest. Meanwhile, distance to faults values remain negative and range from −0.0014 to −0.000016, as described in Figure 5c. Figure 5d represents the variation of slope angle values; negative values were only found in the northwestern part of the study area and range from −0.034 to −0.016, whereas other parts of Penang Island show positive values ranging from 0.00081 to 0.054. Figure 5e shows that the curvature values are constantly negative and range from −0.72 to −0.040, meanwhile slope aspect values range from −0.00033 to 0.003, as described in Figure 5f; slope aspect values are only negative in the southwest. For annual rainfall average, the estimated parameter is not always positive, as values range from −0.0016 to −0.00031 and are scattered in the center of Penang Island and towards the northwest (Figure 5g). Figure 5h depicts the soil series factor, which constantly shows negative values ranging from −2.32 to −0.41. Similarly, the lithology factor in Figure 5i shows that the majority of lithology types deliver negative values ranging from −0.98 to 0.12. Lastly, Figure 5j represents the estimated parameter for land use; values range from −0.29 to 0.073 and the majority of them are negative.

These values represent the significance of each landslide factor across the study area. However, in GWLR, estimate values are not the only indicator of the landslide factors’ significance, and as a result, they should always be associated with the t-values.

Therefore, Figure 6 demonstrates a compilation of 10 integrated maps of the landslide factors considering each t-value. Each variable is distinguished by a 5-color scheme to visualize the significance level in different areas. Figure 6a,c,e,h shows that distance to roads, distance to faults, curvature, and soil series values are negative across the study area, indicating that roads, faults, curvature, and soil series factors have a negative influence on landslide occurrences. Apart from that, negative values for distance to rivers and lithology were scattered across the study area, except for in the southwest and in the northwest, respectively (Figure 6b,i). Meanwhile, Figure 6g describes the annual rainfall average factor, showing that the t-values do not always stay positive, as negative values are found in the center and towards the northwestern part of Penang Island.

Figure 5 and Figure 6 both suggest that both the slope angle and the slope aspect are positively correlated with landslides, as the values increased, thus increasing the probability of landslide occurrences. This indirectly shows that both the slope angle and the slope aspect are the most significant factors influencing landslide occurrences in Penang Island, necessitating comprehensive mitigation actions. Because Penang Island is made up of steep topography above 60 m (exceeding 50 percent of the land area), half of Penang’s land is vulnerable to landslides. As the slope angle increases, shear stress in soil or other unconsolidated material generally increases as well [45]. Even though the slope angle is one of the factors deemed to be significant, it does not always remain a positive influence. The results show that the northwestern part conveys small negative values, indicating that the slope angle is less likely to contribute to landslide occurrences there. Meanwhile, according to [46,47,48], apart from the East-facing slopes, a slope that does not face the sun has a low temperature with a high soil moisture that when it is disturbed, it will tend to influence the landslide occurrences. Therefore, it is necessary to reinforce the slope treatment in other parts of Penang Island, instead of the northwestern part.

5. Conclusions

Several factors of landslides, namely distance to roads, distance to rivers, distance to faults, slope angle, slope aspect, curvature, rainfall annual average, lithology, soil series, and land use have been considered in this study. Based on the results obtained from this study, the landslides in Penang Island are spatially influenced by both the slope angle and the slope aspect. It was mentioned before that half of Penang Island is vulnerable to landslides due to its topography, and thus for urbanization activities, cutting slopes is inevitable. Moreover, a higher inclination of the slope can contribute to a higher gravity force for pulling materials down the slope, thereby increasing the risk of landslides. Therefore, a proper perspective and a thorough understanding of the certain slope conditions have to be established to avoid more landslide occurrences in the future and to ensure the conservation, restoration, and sustainable use of terrestrial and inland freshwater ecosystems, as stated in the 15th SDG. In the meantime, the maps generated can indirectly be useful to authorities in terms of coordinating landslide assessment guidelines, as mentioned by Ministry of Housing and Local Government of Malaysia, as well as by decision makers for selecting suitable locations for future development. We suggest that a future study should highlight an in-depth investigation of the potential use of GWLR at various scales so that the method can be fully utilized in other regions to tackle natural disaster impacts.