Detecting spatially non-stationary and scale-dependent relationships between urban landscape fragmentation and related factors using Geographically Weighted Regression

Applied Geography 31 (2011) 292e302

Contents lists available at ScienceDirect

Applied Geography journal homepage: www.elsevier.com/locate/apgeog

Detecting spatially non-stationary and scale-dependent relationships between urban landscape fragmentation and related factors using Geographically Weighted Regression Jiangbo Gao, Shuangcheng Li* College of Urban and Environmental Sciences, Peking University, Laboratory for Earth Surface Processes, Ministry of Education, Beijing 100871, China

a b s t r a c t Keywords: Geographically Weighted Regression Landscape fragmentation Spatial non-stationarity Scale-dependence Shenzhen City China

Landscape fragmentation is usually caused by many different anthropogenic influences and landscape elements. Scientifically revealing the spatial relationships between landscape fragmentation and related factors is highly significant for land management and urban planning. The former studies on statistical relationships between landscape fragmentation and related factors were almost global and single-scaled. In fact, landscape fragmentations and their causal factors are usually location-dependent and scaledependent. Therefore, we used geographically Weighted Regression (GWR), with a case study in Shenzhen City, Guangdong Province, China, to examine spatially varying and scale-dependent relationships between effective mesh size, an indicator of landscape fragmentation, and related factors. We employed the distance to main roads as a direct influencing factor, and slope and the distance to district centers as indirect influencing factors, which affect landscape fragmentation through their impacts on land use and urbanization, respectively. The results show that these relationships are spatially nonstationary and scale-dependent, indicated by clear spatial patterns of parameter estimates obtained from GWR models, and the curves with a characteristic scale of 12 km for three explanatory variables, respectively. Moreover, GWR models have better model performance than OLS models with the same independent variable, as is indicated by lower AICc values, higher Adjusted R2 values from GWR and the reduction of the spatial autocorrelation of residuals. GWR models can reveal detailed site information on the different roles of related factors in different parts of the study area. Therefore, this finding can provide a scientific basis for policy-making to mitigate the negative effects of landscape fragmentation. Ó 2010 Elsevier Ltd. All rights reserved.

Introduction Landscape fragmentation due to road construction, urbanization, land use/land cover change (LUCC) and other anthropogenic factors leads to more and smaller habitat patches, increased isolation among habitat patches, decreased complexity of patch shape, and higher proportions of edge habitat (Saunders, Mislivets, Chen, & Cleland, 2002). It is considered as one of the most serious threats to the conservation of natural ecosystems as well as the health of agricultural and urbanized ecosystems (Zeng & Wu, 2005). Therefore, the phenomenon of landscape fragmentation has recently attracted more and more attention due to growing applications in biological conservation and ecosystem management (Velázquez, Bocco, & Romero, 2003; Velázquez et al., 2009). However, most studies on landscape fragmentation focus on the depiction of

* Corresponding author. Tel.: þ86 10 62767428; fax: þ86 10 62751187. E-mail addresses: [email protected] (J. Gao), [email protected] (S. Li). 0143-6228/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.apgeog.2010.06.003

landscape pattern characteristics (Feranec, Jaffrain, Soukup, & Hazeu, 2010; Long, Nelson, & Wulder, 2010) and the indication of its impacts on wildlife habitats (Liu, Li, & Li, 2007), biodiversity (Trombulak & Frissell, 2000), ecological processes and functions (Givertz, Thorne, Berry, & Jaeger, 2008), but few of these studies indicate the causes of landscape fragmentation (Munroe, Croissant, & York, 2005), which is highly significant for land management and urban planning. Landscape fragmentation is usually caused by many different human activities such as urbanization and LUCC, and landscape elements such as roads, railways and rivers, whereas a merely qualitative description of the causes is not convincing. Therefore, the quantitative indication of the relationships between landscape fragmentation and its impact factors should be strengthened. In order to effectively characterize the process of landscape fragmentation and its impact factors, the measure of landscape fragmentation should be done firstly. Many landscape fragmentation metrics were proposed to quantify spatial segregation (Li, Chang, Peng, & Wang, 2009). One such metric used in this study

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

is the effective mesh size (meff), which explicitly incorporates the ecological process of animal dispersal into its definition. Effective mesh size is an expression of the probability that any two randomly chosen locations in the landscape are connected, i.e. not separated by barriers such as transportation infrastructure or urban areas (Jaeger, 2000). Spatial data exhibit two properties, i.e. spatial autocorrelation and non-stationarity, which makes it difficult to meet the assumptions and requirements of conventional regression techniques such as Ordinary Least Squares (OLS). Traditional statistical methods can only produce “average” and “global” parameter estimates (Bacha, 2003; Batisani & Yarnal, 2009; Geri, Amici, & Rocchini, 2010), and thus they are unable to deal with spatial autocorrelation existing in the variables. In recent years, a relatively simple, but effective, new technique for exploring spatially varying relationships, called Geographically Weighted Regression (GWR), has been developed (Brunsdon, McClatchey, & Unwin, 2001; Fotheringham, Charlton, & Brunsdon, 2001). GWR allows different relationships to exist at different points in the study area and improves the modeling performance by reducing spatial autocorrelations. In addition, these relationships also greatly depend on scale, which is inherent in natural and man-made processes and patterns (Lü & Fu, 2001). Therefore, local rather than global parameters can be estimated, and spatial non-stationarity can be detected at multi-scales by changing bandwidth of GWR. The objective of this paper, with a case study in Shenzhen City, Guangdong Province, China, was to investigate the applicability of GWR in modeling the relationships between effective mesh size and related factors, and then examine their spatial non-stationarity and the scale-dependence. In this study, effective mesh size calculated using FRAGSTATS software is used as dependent variable, and three impact factors of urban landscape fragmentation are employed as explanatory variables in the regression. Study area We used Shenzhen City as our case study, which lies in the south of Guangdong Province in southern China (Fig. 1), at 22 260 e22 510 N and 113 450 e114 370 E, and is the passageway from mainland China to Hong Kong Special Administrative Region. It has a total terrestrial area of 1952.84 km2 (including all islands) with a NortheSouth span longer than its EasteWest. Vegetation covers 50e80% of the land area. Its topography is relatively higher in the southeast part and lower in the northwest part with elevation ranging from 1 m (coastal zone) to 935 m (Fig. 1). The primary landscape type is hill, with dotted patches of plain. With a southern subtropical monsoon climate, the average annual rainfall is 1938 mm, and the mean annual temperature is 22.4  C. The main soil types are yellow soil, red soil and lateritic red soil according to soil horizontal zonality. Shenzhen Municipality was established in 1979 and was designated a Special Economic Zone in 1980 to attract domestic and foreign investment with tax relief, favorable land development policies, and a skillful and cheap labor force. For only 30 years, Shenzhen City has conglomerated nearby small rural towns and numerous fishing villages with less than 20,000 people into one of the largest cities in China with a population of 8.46 million. Although the opening up of other cities and regions for investment in China in the late 1980s has reduced its comparative advantage, the city has still become one of the most vigorous economies with a gross domestic product (GDP) of 7807 billion RMB in 2008, because of the connections with foreign investors that have been established for three decades (Yeh & Xu, 1996), and the geographical proximity to Hong Kong, which is unrivaled by other Chinese coastal cities (Sui & Zeng, 2001). Owing to the rapid

293

industrialization and urbanization, large areas of natural ecosystems have been converted into construction land use types since 1979, resulting in severe urban landscape fragmentation. Materials and methods The pattern, process and spatial relationships are fundamental issues in geography (Li & Cai, 2005). Landscape patterns interact intensively with ecological processes (Gustafson, 1998) and are usually considered as the results of various ecological processes at multi-scales. Based on sufficient consideration of ecological processes and their key influencing factors that can have important impacts on landscape fragmentation, this paper is intended to indicate spatially varying and scale-dependent relationships between landscape fragmentation and its related factors in Shenzhen City using GWR. Dependent variable: effective mesh size Effective mesh size quantifies landscape fragmentation based on the probability that any two randomly chosen points in the landscape are located in the same non-fragmented patch. It can also be interpreted as the average size of the area that an animal placed randomly in the landscape will be able to access without crossing barriers. The more barriers in the landscape, the lower probability that two points will be connected, and the lower the effective mesh size (Givertz et al., 2008; Jaeger, 2000). Compared with previous landscape metrics such as edge density, which is a function of the length of edge at landscape level regardless of any patch affiliation, effective mesh size is more comprehensive due to the incorporation of ecological process, landscape composition, landscape spatial configuration, and landscape shape into its definition and computation. The spatial pattern analysis program FRAGSTATS (McGarigal & Marks, 1995), developed by Dr. Kevin McGarigal with programming by Eduard Ene and additional programming assistance by Chris Holme, was used to quantify the landscape fragmentation, because it facilitates the spatial analysis of landscape patches and the modeling of the associated attributes. In this study, Landsat Enhanced Thematic Mapper Plus (Landsat/ETMþ) images of Shenzhen (path/row 121-44 and 122-44), acquired on 2nd January 2000, were rectified to match other georeferenced data layers. They were then classified to two different land use types (i.e. build environment and greenspace) (Fig. 2) by supervised classification in ERDAS 8.5. The results of land classification were used to further calculate effective mesh size of land use in FRAGSTATS with a moving window of 60 m radius. The formula is as follows: n X Aij mefff ðjÞ ¼ Aj $ Aj i¼1

!2 ¼

n 1X A2 Aj i ¼ 1 ij

(1)

where n equals the number of unfragmented patches in landscape j, Aij is the size of patch i within landscape j, and Aj is the total area of landscape j. The values of meff range from the cell size when the landscape is maximally subdivided (i.e. every cell is a separate patch), to total landscape area when the landscape consists of a single patch. Fig. 3 is the spatial pattern of effective mesh size retrieved with above-mentioned algorithm and data. Explanatory variables Distance to main roads Road construction, which was intensified in the past 30 years to improve the transportation accessibility in Shenzhen City, has important impacts on ecosystem and landscape fragmentation. Roads directly caused landscape fragmentation as one of several

294

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

Fig. 1. Location of study area. Shenzhen is located in Guangdong Province at the mouth of the Pearl River Delta, and adjacent to Hong Kong in southern China.

Fig. 2. The results of land classification in ShenzhenCity. The Landsat ETM þ images were originally classified to six different land use types (build environment, forest land, garden land, green land, wetland, and cropland), and then the latter five types were reclassified to one type, greenspace.

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

295

Fig. 3. The spatial pattern of effective mesh size of land use calculated with the results of land classification in FRAGSTATS.

landscape elements. Thus, the distance to main roads (abbreviation: DRDS) could reflect and affect the situation of landscape fragmentation. This distance was consequently selected as a direct influencing factor. Furthermore, it is important to investigate the fragmentation characteristics for road management and its surrounding ecosystem (Li, Xu, Zhou, & Wang, 2004). Distance to district centers The natural land covers in built-up areas vary according to the different distances to city centers, so this distance can usually represent and influence the degree of urbanization, and could also affect landscape fragmentation. Shenzhen City has experienced rapid urbanization and LUCC in the past 30 years. Therefore, the distance to district centers (abbreviation: DCTY) was selected as an explanatory variable, which affects landscape fragmentation through its impacts on ecological processes such as LUCC and urbanization. Slope The observed spatial configuration of land use is the collective result of myriad factors that influence the returns, or benefits, to land use (Munroe & York, 2003). Slope is one of the key factors affecting land use, so it could have important impacts on landscape fragmentation. In the study area, most build environment is located in low slope locations while most greenspaces are found in high slope locations (Figs. 1 and 2). Generally, the greater the human activity, the higher the degree of landscape fragmentation. The pattern of landscape fragmentation is thus usually in accordance with the pattern of slope. Therefore, slope was also employed as a predictor that affects landscape fragmentation through its impacts on human activities such as land use. The Pearson correlations between effective mesh size and explanatory variables calculated using SPSS software version 15.0 (SPSS Inc., Chicago, IL, USA) are shown in Table 1. All the related factors exhibit significant positive correlations with landscape fragmentation. Generally, the closer to the built-up area and road, the more fragmented the landscape is, and the lower the effective

mesh size is. Low slope areas are usually in accordance with more human activities, where landscape fragmentation is usually more serious with lower effective mesh size. Therefore, based on the results of Pearson correlations, these three explanatory variables are suitable predictors of effective mesh size over space. In addition, according to the results of one-sample KolmogoroveSmirnov tests in SPSS 15.0, all dependent and explanatory variables meet the condition of normal distribution, which means they are suitable for regression analysis. Geographically Weighted Regression Ordinary Least Squares regression (OLS) and Geographically Weighted Regression (GWR) were used to investigate the relationships between landscape fragmentation and related factors. Since OLS is well known, we will give only a brief introduction for the theoretical background of the GWR model in the next paragraphs. Moreover, steps for data pre-processing and stationary index calculation are also described in this section. Brief description of theoretical background of GWR Depending on the local context in the study area, the relationships between landscape fragmentation and related factors can be different. Differences of distances to district centers and main roads, and slopes usually cause different degrees of landscape Table 1 Pearson correlations between effective mesh size and explanatory variables. Explanatory variable

N

R

p_value

DCTY DRDS Slope

4990 4990 4990

0.610 0.330 0.591

0.01 0.01 0.01

DCTY and DRDS represent the straight line distance to district centers and main roads, respectively.N is the number of sampling data; R is the Pearson correlation; and p_value is the pseudo-significance level, which is computed using a randomization algorithm.

296

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

fragmentation. OLS, one of global spatial models, which examines the average state, can be misleading when describing phenomena that vary over space (Clement, Orange, Williams, Mulley, & Epprecht, 2009). GWR is an extension of traditional standard regression techniques such as OLS because it allows local rather than global parameter estimates (Fotheringham et al., 2001). These estimates, showing how a relationship varies over space, can help to examine the spatial pattern of the local estimates to get some understanding of hidden possible causes of this pattern (Fotheringham, Brunsdon, & Charlton, 2002). Furthermore, the GWR model is easy to use because it was added to the ArcToolbox in ArcGIS 9.3 (ESRI Inc., 1999e2008). The conventional global regression can be expressed as:

b y i ¼ b0 þ

X

bk xik þ 3i ;

(2) Results

k

where b y i is the estimated value of the dependent variable at location i, b0 represents the intercept, bk expresses the slope coefficient for independent variable xk, xik is the value of the variable xk at location i, and 3i denotes the random error term for location i. In this equation, the estimates of the model parameters are assumed to be spatially stationary. The GWR model extends conventional global regression by generating a local regression equation for each observation, and the above model can be rewritten as:

b y i ¼ b0 ðmi ; ni Þ þ

X

bk ðmi ; ni Þxik þ 3i ;

(3)

k

where (mi, ni) denotes the coordinate location of the ith point, b0(mi, ni) is the intercept for location i, bk(mi, ni) represents the local parameter estimate for independent variable xk at location i. Parameter estimates in GWR are obtained by weighting all observations around a specific point i based on their spatial proximity to it. The observations closer to point i have higher impact on the local parameter estimates for the location, and are weighted more than data far away. The parameters are estimated from:

b b ðm; nÞ ¼



X T Wðm; nÞX

1

X T Wðm; nÞy;

(4)

b (m, n) represents the unbiased estimate of b, W(m, n) is the where b weighting matrix which acts to ensure that observations near to the specific point have bigger weight value. The weighting function, called the kernel function, can be stated using the exponential distance decay form: wij ¼ exp

d2ij b2

The measurement of spatial non-stationarity In this paper, we calculated the stationary index using Matlab 7.4 (The MathWorks Inc. 1984e2007), which was proposed by Brunsdon, Fotheringhm, and Charlton (1998) and Osborne, Foody, and Suárez-Seoane (2007) (Brunsdon et al., 1998; Osborne et al., 2007). This index is designed to measure the spatial nonstationarity, and values smaller than one indicates stationarity (Charlton, Fotheringhm, & Brunsdon, 2003). The calculation included three steps: Firstly, we computed the interquartile range of standard error for GWR coefficients for each explanatory variable using Matlab function iqr; secondly, twice standard error of the global regression coefficient was obtained with Matlab function glmfit; finally, the ratio of these two factors was used as the stationary index.

! (5)

where wij represents the weight of observation j for location i, dij expresses the Euclidean distance between points i and j, and b is the kernel bandwidth. If observation j coincides with i, the weight value is one. If the distance is greater than the kernel bandwidth, the weight will be set to zero. Data pre-processing To meet the requirement of ArcGIS-based GWR, all GIS layers for dependent and independent variables, which were raster data at 30m spatial resolution initially, were converted to vector formats. Three steps of data conversion are as follows: (1) 5000 points in the study area were selected randomly by using Hawth’s Analysis Tools for ArcGIS (http://www.spatialecology.com/htool/). (2) Values were extracted from raster data layers to 5000 points using ArcToolbox in ArcGIS9.3 (ESRI Inc., 1999e2008). (3) 4990 points were used to GWR model after removing 10 NoData cases.

Although it is clear that human activities influence landscape fragmentation, the relationships between landscape fragmentation and related influencing factors and their spatial heterogeneity often remain unclear. Furthermore, these relationships are usually sensitive to scale due to the mismatch between observation scale and intrinsic scale. Observation scale is the measurement of spatial pattern and processes being studied, while intrinsic scale is the inherent organization of the natural phenomena (Chen, Liu, Lü, Feng, & Fu, 2008). Only when the observation scale is in accordance with the intrinsic scale can these relationships be exactly indicated (Wu, 2004). Scale-dependence of spatial relationships Stationary index and two diagnostics (Residualsqaures and Sigma) of GWR models at multi-scales for three explanatory variables are calculated. Residualsqaures is the sum of the squared residuals in the model, while the residual is the difference between an observed value and its estimated value returned by the GWR model. Sigma is the square root of the normalized residual sum of squares, where the residual sum of squares is divided by the effective degrees of freedom of the residual. The smaller the values of these two statistics, the closer the fit of the GWR model to the observed data. As shown in Fig. 4, the stationary index declines rapidly with a coarsening of the scale, especially for the variable DRDS, and the declines of the stationary index for all predictors are largest at a spatial scale of 4 km, which means that the dominant operating scales of each independent variable on effective mesh size are located within 4 km. Compared with the other two variables, DRDS has not only the bigger decline gradient, but also a relative larger stationary index, which indicates that the distance to main roads influences the effective mesh size at broader spatial scale range. Although the decline slopes of curves become quite flat at a spatial scale of above 12 km, only the variable Slope became stationary within the spatial scale range of 15 km, whose stationary index is smaller than one at a spatial scale of 4.3 km. Likewise, Koenker (BP) Statistic (Koenker’s studentized BrueschePagan statistic) also indicated statistically significant heteroscedasticity and/or nonstationarity at the level of p < 0.05 for the other two variables, DRDS and DCTY, according to the statistical result of OLS model. The results of two diagnostics at multi-scales, displaying nearly opposite variation trends compared to Fig. 4, show that all the rising slopes become quite flat at a spatial scale of above 12 km. The variable DCTY has a better fit than other explanatory variables within the spatial scale range. It can be concluded that the interrelations between landscape fragmentation and related factors exhibit significant spatial scale-

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

297

Fig. 4. Stationary indexes at multi-scales for three explanatory variables. The stationary index is a ratio between the interquartile range of standard error for GWR coefficients and twice standard error from global regression analysis such as OLS, and values smaller than one indicates stationary.

dependence. The spatial scale of 12 km, the characteristic scale of all curves, above which curves of the stationary index and two diagnostics for predictors become quite flat, can be regarded as the intrinsic scale of these relationships. This scale can be used as the operating scale in this study (i.e. the bandwidth for GWR). Thus, GWR with a bandwidth of 12 km is the suitable technique for modeling the relationship between effective mesh size and related factors in the study area, rather than OLS with global constant regression coefficients. Comparisons between OLS and GWR models Comparison of model performance between OLS and GWR In this study, we compared AICc and Adjusted R2 values from GWR with those from OLS to investigate whether GWR models have better model performance than the corresponding OLS models. AICc (corrected Akaike’s Information Criterion) is a measurement of model performance which is useful for comparing different regression models. The model with lower AICc value means a better fit to the observed data and better model performance. Usually, a decrease of AICc values lower than three indicates the model with the lower AICc is held to be better. R-Squared (R2) measures goodness of fit, whose value varies from 0 to 1. Adjusted R2 improves R2 by normalizing the numerator and

denominator on the basis of their degrees of freedom to give more real model fits, and it is almost always smaller than the R2 value. Higher Adjusted R2 indicates that the regression model can account for more dependent variable variance. The AICc and Adjusted R2 of GWR and OLS models are shown in Table 2. The tabular data from second row to seventh row in Table 2 are the results of GWR and OLS models between effective mesh size and each explanatory variable. The data of last two rows are the results of GWR and OLS models with all combined explanatory variables. In all cases, the GWR models outperform the OLS models, Table 2 Comparison of two diagnostics (AICc and Adjusted R2) between OLs and GWR models. Adjusted R2

Explanatory variables

AICc

DCTY

AICO AICG

9051 7729

Adjusted R20 Adjusted R2G

0.372 0.519

DRDS

AICO AICG

10800 9282

Adjusted R20 Adjusted R2G

0.109 0.344

Slope

AICO AICG

9233 8642

Adjusted R20 Adjusted R2G

0.349 0.423

All variables

AICO AICG

8105 7062

Adjusted R20 Adjusted R2G

0.481 0.581

AICO is the AICc for OLS model; AICG is the AICc for GWR model.

298

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

Table 3 Comparison of Moran’s I of residuals between OLS and GWR. Explanatory variables

IO

p_value

IG

p_value

DCTY DRDS Slope All variables

0.49 0.55 0.39 0.39

0.01 0.01 0.01 0.01

0.36 0.39 0.34 0.28

0.01 0.01 0.01 0.01

IO is the Moran’s I for OLS model; IG is the Moran’s I for GWR model; and p_value is the pseudo-significance level.

characterized by the higher Adjusted R2 values and the lower AICc values of the GWR models. The local models echo and amplify the relationships observed by the global models, so the performance of the local models represents a complex mix of the model’s overall performance and the relative degree of non-stationarity (Lafary, Gatrell, & Jensen, 2008). Therefore, the relative performance of the global and local models illustrate that the parameters display non-stationarity, and that these issues can be better understood when the statistical analyses are ‘scaled down’ enabling case-wise examination of the parameter diagnostics. Comparison of spatial autocorrelations of residuals between OLS and GWR Autocorrelation is an ecological data expression for the lack of independence between pairs of observations at given distances apart, both in time and space (Legendre, 1993). Generally, the spatial randomness of residuals and the cumulative residual amount are employed as the diagnostic measures of the regression model. Statistically significant clustering of high and/or low residuals indicates that at least one key explanatory variable, which could effectively capture the inherent spatial structure in the dependent variable, is missing from the model. In order to compare the ability to deal with spatial autocorrelation between OLS and GWR, global and local Moran indexes (Moran’s I) of residuals from each of the OLS and GWR models were computed using the ArcToolbox in ArcGIS 9.3 (ESRI Inc., 1999e2008) and the GSþ 7.0

version (Gamma Design Software, 2004), respectively. Moran’s I is commonly used as an indicator of spatial autocorrelation, and its values range from 1 to 1. The larger the absolute value of Moran’s I, the more significant the spatial autocorrelation. A value of 0 means perfect spatial randomness. The results of global Moran’s I statistics on the residuals from OLS and GWR models are shown in Table 3. The tabular data from the second row to the fourth row in Table 3 are the results of GWR and OLS models between effective mesh size and each explanatory variable. The data of the last row are the results of GWR and OLS models with all combined explanatory variables. Significant positive spatial autocorrelations are found for all the OLS models, characterized by Moran’s I ranging from 0.39 to 0.55 and p < 0.01, and all the GWR models, characterized by Moran’s I ranging from 0.28 to 0.39 and p < 0.01. The results show that the GWR models produced smaller global Moran’s I than OLS models with the same independent variable, indicating that GWR models improve the reliability of the relationships by reducing the spatial autocorrelations in residuals. In order to more elaborately detect the spatial structure of residuals from OLS and GWR models, the local Moran’s I was calculated at multi-scales. The results show that all the local Moran’s I statistics on the residuals from GWR models within the spatial scale range are smaller than those from OLS models. Besides, no significant spatial autocorrelation was found over 2500 m for the GWR model, while spatial autocorrelations of residuals from OLS models are relatively significant. The results indicate that the OLS assumption of residual independence is not met, and the GWR model is an improvement of the OLS model. However, the application of GWR requires a careful use of diagnostic parameters such as Moran’s I statistics. If an OLS model has spatial autocorrelation problem, GWR can help reduce it. On the other hand, if an OLS model does not have this problem, an application of GWR may increase spatial autocorrelation (Tu & Xia, 2008). The spatial autocorrelations of residuals from OLS models might be caused by the spatial distribution of effective mesh size, which is

Fig. 5. Cluster map of effective mesh size. Cluster definitions: HigheHigh: High values surrounded by high values; HigheLow: High values surrounded by low values; LoweHigh: Low values surrounded by high values; LoweLow: Low values surrounded by low values.

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

299

Fig. 6. Spatial variation of regression outputs from the GWR model for effective mesh size and distance to district centers (DCTY): (a) Slope Parameter for DCTY; (b) Local R-squared (R2) for DCTY; (c) Standardized residuals for DCTY.

not adequately accounted for by OLS models. The results of the local indicators of spatial association (LISA) analysis on effective mesh size are illustrated in Fig. 5. We reduced data points to 2000 in order to make patterns more clearer. As shown in Fig. 5, in points with significant cluster, most high-value points are surrounded by highvalue points, and most low-value points are surrounded by lowvalue points. The cluster map of effective mesh size shows several significant spatial clusters and clear spatial distribution patterns of these clusters. Most LoweLow clusters are located in the west of Shenzhen City (the highly urbanized central city area), and most HigheHigh clusters are located in the southeast (relatively lessurbanized rural areas). Moreover, significant positive spatial autocorrelations are also found for effective mesh size by the calculation of global Moran’s I (Moran’s I ¼ 0.50, p < 0.01). Spatial heterogeneity of relationships between landscape fragmentation and related factors Maps of slope parameters (b coefficients), local R2, and standardized residuals (stdresid) obtained from GWR models provide a simple way to detect the spatial varying relationships between effective mesh size and related factors. Local R2, ranging from 0 to 1, indicates how well the local regression model fits the observations, and the local models with low values perform poorly. Residuals are the differences between the observed y values and the predicted y values. Standardized residuals have a mean of zero and a standard deviation of 1. We also computed the fractal dimension of slope parameter, which provides a means of quantitatively expressing the complexity of the data. Generally, high value indicates a complex series (Phillips, 1986).

The maps of slope parameters, local R2, and standardized residuals from GWR model for effective mesh size and DCTY are shown in Fig. 6. DCTY has a significant positive correlation with the effective mesh size (R ¼ 0.610, p < 0.01; Table 1), which indicates that higher effective mesh size values are related to greater distances to district centers. Similarly, significant positive correlations are also observed for the entire study area from the GWR model (Fig. 6a). However, a clear spatial pattern can be identified from the results of the GWR model (Fig. 6a and b). More significant correlations and higher local R2 values are located in the northwest, and less significant correlations and lower local R2 values are found in the southeast. DCTY can explain more than 40% of the spatial variation in effective mesh size in most built-up areas but can explain less than 30% in most greenspace areas (Fig. 6b), while the OLS model can account for only 37% of the variance. This result means that the distance to built-up area has a more important impact on the effective mesh size in highly urbanized areas than it does in less-urbanized areas. Although the positive values of standardized residuals indicating over-prediction appear in greenspace areas and negative values indicating under-prediction appear in build environment areas, the standardized residuals from the GWR model show a relatively random pattern of over-/underprediction (Fig. 6c) compared with the other two maps. As opposed to the significant positive correlation between effective mesh size and DRDS obtained from OLS model (R ¼ 0.330, p < 0.01; Table 1), both significant positive and negative correlations are observed by the GWR model (Fig. 7a). Significant spatial non-stationarity can be also observed in the maps of parameter estimates and local R2 (Fig. 7a and b). Stronger positive relationships between effective mesh size and the distance to main roads

Fig. 7. Spatial variation of regression outputs from GWR model for effective mesh size and distance to main roads (DRDS): (a) Slope Parameter for DRDS; (b) Local R-squared (R2) for DRDS; (c) Standardized residuals for DRDS.

300

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

Fig. 8. Spatial variation of regression outputs from the GWR model for effective mesh size and Slope: (a) Slope Parameter for Slope; (b) Local R-squared (R2) for Slope; (c) Standardized residuals for Slope.

and higher local R2 are located in central part of the study area, while negative and weaker relationships and local R2 are found in the west and southeast. The spatial pattern of standardized residuals from the GWR model for effective mesh size and DRDS (Fig. 7c) is highly similar to Fig. 6c. Fig. 8 shows the maps of slope parameters, local R2, and standardized residuals from the GWR model for effective mesh size and slope. Similar to the result of the OLS model for the correlation between effective mesh size and slope (R ¼ 0.591, p < 0.01; Table 1), significant positive correlations are also found for all the sampling sites in the GWR model (Fig. 8a). However, both Fig. 8a and b represent clear clusters. More significant correlations are located in the west and northeast of the study area while lower correlations are found in the southeast (Fig. 8a). Most of the higher local R2 values can be found in the southern part (Fig. 8b). In addition, the spatial pattern of standardized residuals (Fig. 8c) is also highly similar to Fig. 6c. Moreover, the fractal dimensions of slope parameters in Figs. 6 and 7, and Fig. 8 show that the series of coefficients for variable DCTY (fractal dimension ¼ 1.654, R2 ¼ 0.872) are the most complex compared with DRDS (fractal dimension ¼ 1.519, R2 ¼ 0.880) and Slope (fractal dimension ¼ 1.305, R2 ¼ 0.960). Conclusions and discussion Conclusions Scientifically revealing the spatial relationships between landscape fragmentation and related factors is highly significant for land management and urban planning. Shenzhen has experienced the rapid expansion of desakota regions in the past 30 years, and has completed the process of urbanization in 2005 due to rapid industrialization, which not only dramatically increased the builtup areas, but also caused landscape fragmentation in numerous ways. Thus, in this study, we detected spatial relationships between landscape fragmentation and related factors using Geographically Weighted Regression (GWR) with a case study in Shenzhen City, Guangdong Province, China. The most interesting finding of this study is that the spatially non-stationary relationships between landscape fragmentation and related factors are scale-dependent in the study area, indicated by the curves with the characteristic scale of 12 km for three explanatory variables. The results also show that only the variable Slope became stationary within the spatial scale range of 15 km, and the variable DCTY has a better fit than other explanatory variables within the spatial scale range. Our results show that GWR models represent a significant improvement of model performance over OLS models with the same predictor, as is indicated by lower AICc values, higher

Adjusted R2 values from GWR and the reduction of the spatial autocorrelation of residuals. Moreover, it can be found that Moran’s I statistics on the residuals at multi-scales from GWR rise rapidly with a coarsening of the scale. Therefore, the reduction of correlations may be caused by the regional paradigm of the GWR model, which will be studied in future. Clear spatial patterns of regression parameters can be also identified from the results of the GWR models. The former studies on the statistical relationship between landscape fragmentation and related factors were almost global, which deviated from real physical process and pattern. In this study, OLS models show that all explanatory variables have significant positive correlations with landscape fragmentation, the results from GWR models, however, show that their relationships are not constant over space (and are even negative in some locations). In addition, maps of local R2 and standardized residuals obtained from GWR models also represent spatial heterogeneity. Discussion Two basic assumptions of global regression models such as OLS are that model residuals have constant variance (homoscedasticity), and are uncorrelated with each other (no autocorrelation). Violations of the assumptions will reduce the efficiency of the regression and mislead the interpretation of model results (Hamilton, 1992). Actually, landscape fragmentation and its causal factors are usually location-dependent and autocorrelated, because many natural processes exhibit highly spatial heterogeneity and autocorrelation, which make it difficult to meet the assumptions and requirements of the OLS model. Our study detected the spatial relationships between landscape fragmentations and related influencing factors using both OLS and GWR models. The results show that GWR can serve as a new approach for effectively exploring spatially varying relationships and can reduce spatial autocorrelation. GWR models revealed the detailed site information on the different roles of related factors in different parts of the study area, which improves the model ability to explain the local situation of landscape fragmentation, and is helpful for developing more effective plans and policies to mitigate the negative effects of landscape fragmentation. In order to avoid the potential multicollinearity among explanatory variables, at first, only one explanatory variable was used for OLS and GWR models to analyze its association with landscape fragmentation. And then all explanatory variables were combined for models to compare the overall performance of goodness-of-fit with that of single independent variable models. Condition number can evaluate local collinearity. Results are unstable in the presence of local collinearity as indicated by

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

a condition number greater than 30. In this study, all the condition numbers from GWR with all explanatory variables, ranging from 3.24 to 5.32, are smaller than 30, which means that there is no strong local collinearity, and the results are reliable. Although spatially varying coefficient technique is prone to multicollinearity of local coefficients so that the coefficients could be correlated even when there is not collinearity among variables (Griffith, 2008), most researchers concur that GWR can be reliably used as an exploratory technique to understand how models may function differently across regions (Ogneva-Himmelberger, Pearsall, & Rakshit, 2009). Although the effective mesh size is a good indicator to quantify landscape fragmentation, patch-based landscape metrics are usually biased by the boundaries and the extent of a reporting unit (moving windows with 60m radius in this study) if the boundaries segregate patches. Fortunately, a modified effective mesh size calculation accounting for cross-boundary connections has been presented in recent advances of landscape metric theory (Moser, Jaeger, Tappeiner, Tasser, & Eiselt, 2007). Landscape pattern is also sensitive to scale, so it is necessary to correlate the observation scale precisely with the intrinsic scale to address the depiction of landscape pattern exactly. Thus, in the future we will study the scale-dependence of the modified effective mesh size to further refine the pattern analysis in landscape ecology and the regression models. In GWR models, all variables are assumed to be spatially nonstationary, which indicates that relationships between the dependent variable and the explanatory variables vary over space. However, in some cases, portions of regression coefficients are constant, while other coefficients are functions of geographical locations. One effective approach is to extend GWR to mixed geographically weighted regression (MGWR) (Mei, Wang, & Zhang, 2006), but the issue of how to distinguish global variables from local variables usually makes it difficult to operate. We believe that the calculations of stationary indexes at multi-scales can provide a simple and effective way to settle this issue, which will be further discussed in the next paper. In the future, we will extent one-dimensional analysis to multidimensional analysis to identify the dynamic causes of landscape fragmentation while considering both composition and configuration of the landscape components, which are measured by the number and distribution of patches or distinct areas of the same land cover type. We will also add more ecological processes to GWR models to provide a more comprehensive basis for ecosystem management and biodiversity conservation. Acknowledgements The research is supported by the National Natural Science Foundation of China (40635028 and 40771001). References Bacha, C. J. C. (2003). The determinants of reforestation in Brazil. Applied Economics, 35, 631e639. Batisani, N., & Yarnal, B. (2009). Urban expansion in Centre County, Pennsylvania: spatial dynamics and landscape transformations. Applied Geography, 29, 235e249. Brunsdon, C., Fotheringhm, A. S., & Charlton, M. (1998). Geographically weighted regression-modeling spatial non-stationarity. Statistician, 47, 431e443. Brunsdon, C., McClatchey, J., & Unwin, D. J. (2001). Spatial variations in the average rainfall-altitude relationship in Great Britain: an approach using geographically weighted regression. International Journal of Climatology, 21, 455e466. Charlton, M., Fotheringhm, A. S., & Brunsdon, C. (2003). GWR 3: Software for geographically weighted regression. England: Spatial Analysis Research Group, Department of Geography, University of Newcastle upon Tyne. Chen, L. D., Liu, Y., Lü, Y. H., Feng, X. M., & Fu, B. J. (2008). Pattern analysis in landscape ecology: progress, challenges and outlook. Acta Ecologica Sinica, 28 (11), 5521e5531, (in Chinese).

301

Clement, F., Orange, D., Williams, M., Mulley, C., & Epprecht, M. (2009). Drivers of afforestation in Northern Vietnam: assessing local variations using geographically weighted regression. Applied Geography, 29, 561e576. Feranec, J., Jaffrain, G., Soukup, T., & Hazeu, G. (2010). Determining changes and flows in European landscapes 1990e2000 using CORINE land cover data. Applied Geography, 30, 19e35. Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically weighted regression: The analysis of spatially varying relationships. West Sussex: John Wiley & Sons. Fotheringham, A. S., Charlton, M. E., & Brunsdon, C. (2001). Spatial variations in school performance: a local analysis using geographically weighted regression. Geographical and Environmental Modelling, 5, 43e66. Geri, F., Amici, V., & Rocchini, D. (2010). Human activity impact on the heterogeneity of a Mediterranean landscape. Applied Geography, 30, 370e379. Givertz, E. H., Thorne, J. H., Berry, A. M., & Jaeger, J. A. G. (2008). Integration of landscape fragmentation analysis into regional planning: a statewide multiscale case study from California, USA. Landscape and Urban Planning, 86, 205e218. Griffith, D. A. (2008). Spatial-filtering-based contribution to a critique of geographically weighted regression (GWR). Environmental and Planning A, 40 (11), 2751e2769. Gustafson, E. J. (1998). Quantifying landscape spatial pattern: what is the state of the art? Ecosystem, 1, 143e156. Hamilton, L. C. (1992). Regression with graphics: A second course in applied statistics. Belmont, California: Duxbury Press. Jaeger, J. A. G. (2000). Landscape division, splitting index, and effective mesh size: new measures of landscape fragmentation. Landscape Ecology, 15, 115e130. Lafary, E. W., Gatrell, J. D., & Jensen, R. R. (2008). People, pixels and weights in Vanderburgh County, Indiana: toward a new urban geography of humaneenvironment interactions. Geocarto International, 23(1), 53e66. Legendre, P. (1993). Spatial autocorrelation: trouble or new paradigm? Ecology, 74 (6), 1659e1673. Li, S. C., & Cai, Y. L. (2005). Some scaling issues of geography. Geographical Research, 24(1), 11e18, (in Chinese). Li, S. C., Chang, Q., Peng, J., & Wang, Y. L. (2009). Indicating landscape fragmentation using LeZ complexity. Ecological Indicators, 9, 780e790. Li, S. C., Xu, Y. Q., Zhou, Q. F., & Wang, L. (2004). Statistical analysis on the relationship between road network and ecosystem fragmentation in China. Progress in Geography, 23(5), 78e85, (in Chinese). Liu, H. Y., Li, Z. F., & Li, X. M. (2007). Effects of wetland landscape fragmentation on habitats of oriental white storks-a case study on northeastern Sanjiang plain, China. Journal of Natural Resources, 22(5), 817e823, (in Chinese). Long, J. A., Nelson, T. A., & Wulder, M. A. (2010). Characterizing forest fragmentation: distinguishing change in composition from configuration. Applied Geography, 30, 426e435. Lü, Y. H., & Fu, B. J. (2001). Ecological scale and scaling. Acta Ecologica Sinica, 21(12), 2096e2105, (in Chinese). McGarigal, K., & Marks, B. J. (1995). FRAGSTATS: Spatial pattern analysis program for quantifying landscape structure. USDA For. Serv. Gen. Tech. Rep. PN-351. Mei, C. L., Wang, N., & Zhang, W. X. (2006). Testing the importance of the explanatory variables in a mixed geographically weighted regression model. Environmental Planning A, 38, 587e598. Moser, B., Jaeger, J. A. G., Tappeiner, U., Tasser, E., & Eiselt, B. (2007). Modification of the effective mesh size for measuring landscape fragmentation to solve the boundary problem. Landscape Ecology, 22, 447e459. Munroe, D. K., Croissant, C., & York, A. M. (2005). Land use policy and landscape fragmentation in an urbanizing region: assessing the impact of zoning. Applied Geography, 25, 121e141. Munroe, D. K., & York, A. M. (2003). Jobs houses, and trees: changing regional structure, local and land-use patterns, and forest cover in southern Indiana. Growth and Change, 34(3), 299e320. Ogneva-Himmelberger, Y., Pearsall, H., & Rakshit, R. (2009). Concrete evidence & geographically weighted regression: a regional analysis of wealth and the land cover in Massachusetts. Applied Geography, 29, 478e487. Osborne, P. E., Foody, G. M., & Suárez-Seoane, S. (2007). Non-stationarity and local approaches to modeling the distributions of wildlife. Diversity and Distributions, 13, 313e323. Phillips, J. D. (1986). Measuring complexity of environmental gradients. Plant Ecology, 64(22), 95e102. Saunders, S. C., Mislivets, M. R., Chen, J., & Cleland, D. T. (2002). Effects of roads on landscape structure within nested ecological units of the Northern Great Lakes Region, USA. Biological Conservation, 103, 209e225. Sui, D. Z., & Zeng, H. (2001). Modeling the dynamics of landscape structure in Asia’s emerging desakota regions: a case study in Shenzhen. Landscape and Urban Planning, 53, 37e52. Trombulak, S. C., & Frissell, C. A. (2000). Review of ecological effects of roads on terrestrial and aquatic communities. Conservation Biology, 14, 18e30. Tu, J., & Xia, Z. G. (2008). Examining spatially varying relationships between land use and water quality using geographically weighted regression I: model design and evaluation. Science of the Total Environment, 407, 358e378. Velázquez, A., Bocco, G., & Romero, F. J. (2003). A landscape perspective on biodiversity conservation: the case of Central Mexico. Mountain Research and Development, 23(3), 240e246.

302

J. Gao, S. Li / Applied Geography 31 (2011) 292e302

Velázquez, A., Cué-Bär, E. M., Larrazábal, A., Sosa, N., Villaseñor, J. L., McCall, M., et al. (2009). Building participatory landscape-based conservation alternatives: a case study of Michoacán, Mexico. Applied Geography, 29, 513e526. Wu, J. G. (2004). Effects of changing scale on landscape pattern analysis: scaling relations. Landscape Ecology, 19, 125e138.

Yeh, A. G. O., & Xu, X. (1996). Globalization and the urban system in China. In F. Lo, & Y. Yeung (Eds.), Emerging World cities in Pacific Asia. Hong Kong: United nations University Press. Zeng, H., & Wu, X. B. (2005). Utilities of edge-based metrics for studying landscape fragmentation. Computers, Environment and Urban Systems, 29, 159e178.