- Email: [email protected]
Geostatistical techniques for incorporating spatial correlation into land use change models
International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446 www.elsevier.com/locate/jag
Geostatistical techniques for incorporating spatial correlation into land use change models Ademola K. Braimoh a,b,*, Takashi Onishi b a
United Nations University Institute of Advanced Studies, International Organizations Center, Yokohama 220-8502, Japan b Research Center for Advanced Science and Technology, University of Tokyo, 4-6-1 Komaba, 153-8904, Japan Received 21 March 2006; accepted 19 February 2007
Abstract Land use modeling requires large amounts of data that are typically spatially correlated. This study applies two geostatistical techniques to account for spatial correlation in residential land use change modeling. In the first approach, we combined generalized linear model (GLM) with indicator kriging to estimate the posterior probability of residential development. In the second approach, generalized linear mixed model (GLMM) was used to simultaneously model spatial correlation and regression fixed effects. Spatial agreement between actual and modeled land use change was higher for the GLM incorporating indicator kriging. The GLMM produced more reliable estimates and could be more useful in analyzing the effects of driving factors of land use change for land use planning. # 2007 Elsevier B.V. All rights reserved. Keywords: Generalized linear mixed model; Indicator kriging; Land use change; Lagos; Land use planning
1. Introduction Land use/cover change (LUCC) is one of the most important human impacts affecting ecosystems and environmental processes at local, regional and global scales. Such effects include land degradation and fragmentation, biodiversity loss, alteration in atmospheric composition and climate change (Vitousek et al., 1997). Major changes may yet be expected in the years to come as a result of population growth, and socio-economic and political developments. An understanding of the determinants of LUCC is therefore important for proper land use planning and the
* Corresponding author. Present address: Global Land Project, Sapporo Nodal Office, Hokkaido University, N9W8, Sapporo 060-0809, Japan. Tel.: +81 11 706 4531; fax: +81 11 706 4534. E-mail address: [email protected] (A.K. Braimoh). 0303-2434/$ – see front matter # 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jag.2007.02.005
development of sustainable environmental management policies. Models are required to formalize knowledge about land use change. They are important tools to explore the interactions between land use dynamics and the driving factors. By linking data and theory using formal equations, they help to improve our understanding of the processes underlying land use change. A land use model typically consists of three components: multitemporal land use maps derived from remotely sensed data, a multivariate function of the hypothesized determinants of change, and the resulting prediction map of land use change (Lambin, 1994). Land use models serve three main purposes. First, they help to identify the set of variables that are most important in explaining land use change. Second, they reveal the relative importance of the significant variables. The magnitude of the effect of factors on land change is important because it may indicate areas where policy
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
intervention is most urgently needed, and whether local action to mitigate adverse effect of land change will be sustainable. Last, land use models may be useful in predicting the location of land use change in the immediate future. The ability to predict vulnerable places to land use change helps to mitigate the adverse effects of land change. Land use modeling often requires biophysical and social data that are typically spatially correlated. Spatial correlation is the tendency of objects near in space to be more similar than those further apart. It is a nuisance in conventional statistics, but in land use modeling, it
439
reveals the pattern, and possibly the processes associated with the pattern of land use. Spatial correlation in land use derives from two sources (Overmars et al., 2003). The first source is the distribution of the biophysical and social variables that determines land use pattern and change. The second source is spatial proximity (interaction) between land use types. Spatial interaction in land use might be due to individual choices in land use location, or due to spatial policies in the allocation of land. For instance, the proximity to publicly supplied infrastructure, employment locations or the need for social interaction might
Fig. 1. The study area.
440
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
encourage the clustering of residential land use in an area. Similarly, spatial planning policies may encourage the concentration of business premises in designated commercial centers. Conversely, down-zoning policies may restrict the extent and intensity of urban development in an area. It is therefore essential to account for spatial correlation in land use change models. Ignoring spatial correlation in land use modeling may result in false detection of significant relationships and biased model parameter estimates (Griffith and Layne, 1999). The objectives of this study are to apply and compare two geostatistical techniques for incorporating modeling residential land change. The first method combines logistic regression with Indicator kriging (IK), whereas the second method concurrently models spatial correlation and regression fixed effects. The study applies to residential land use change in Lagos, Nigeria. 2. The study area The study area of approximately 2910 km2 lies between latitudes 68260 and 68500 N and stretches between longitudes 38090 and 38460 E (Fig. 1). Lagos is comprised of depositional landform features: wetlands, barrier islands, beaches, low-lying tidal flats and estuaries. It is characterized by a wet equatorial climate with mean annual rainfall above 1800 mm, mean monthly temperature of 30 8C, and relative humidity ranging from 80 to 100%. The terrain is relatively flat with a mean elevation of about 24 m. The physical nature of the terrain determines the pattern of settlement, and markedly limits the supply of land for urban development. Lagos is the dominant center for non-agricultural production, distribution and business services in Nigeria. Its current population of over 13 million is estimated to grow annually at the rate of 5.6% (UN, 2002). Housing is grossly inadequate to cater for the increasing population. This problem is related to a number of social and institutional factors such as the skewed distribution of private land, the high costs of undeveloped land and housing construction, and weak planning regulations. Thus, new settlements, even if substandard, are only affordable and are springing up at the periphery of the Lagos metropolis. Despite the importance of Lagos in the socio-economic development of Nigeria and for the West African sub-region, there is limited study of the factors responsible for the development of its spatial pattern of land use. Such a dearth of knowledge hinders the development of policies for sustainable land use planning and management of the city.
3. Methods 3.1. Data sources Land use/cover maps of the study area for 1984 and 2000 prepared from Landsat TM images were used. Land cover classification was performed using spectral angle mapping technique (Kruse et al., 1993). Spectral endmembers for the classification were derived from existing land use maps and intensive field survey in May 2005. Fig. 2 shows the maps for four land classes: residential, industrial/commercial, non-urban land, and water; and their relative proportions in 1984 and 2000. The dependent variable for this study, residential development between 1984 and 2000, was extracted from the maps within a geographical information system. Eight spatially explicit independent variables hypothesized to affect residential development were also developed within a GIS (Table 1). These include travel times to major roads and the central business district, Euclidean distances from waterbodies and publicly supplied water stations, population potential in 1984, elevation, frequency of non-urban land within 3 3 neighborhoods in 1984, and number of different land classes within 3 3 neighborhoods in 1984. All variables were mapped at a resolution of 90 m, dictated by the grain size of the digital elevation model. A total
Fig. 2. Land use/cover maps and proportions.
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
441
Table 1 Independent variables list Variable name/description
Abbreviation
Minimum
Mean
Maximum
Elevation (m) Distance from rivers and other waterbodies (m) Distance from waterworks (m) Travel time to Lagos Island, the Central Business District (min) Travel time to major roads (min) Population potential in 1984 Number of non-urban land within 270 m neighborhood Number of different classes within 270 m neighborhood
ELEV WATER WATERWKS CBD
36 0 0 0
24 6,665 13,026 107
115 26,364 36,495 179
ROADS POP84 NURB NDC
0 334 0 1
of 2668 points were randomly selected, out of which 2080 points were used for modeling and the remaining 588 points for validation.
56 6,404 5.99 1.27
23 6532 8337 32 33 8046 3.94 0.47
Coefficient of variation (%) 96 98 64 30 59 126 66 37
maximum likelihood algorithm with equal prior probabilities. After back transformation, the result of the logistic regression may be expressed in terms of conditional probability as
3.2. Generalized linear model (GLM)
ˆ
Generalized Linear Model relaxes the independent and identically distributed assumption commonly used for normal data. It has three distinct parts: a random component which is the qualitative variable being modeled, a systematic or fixed effect component which refers to the independent variables that could be either categorical or continuous; and a link function which is used to relate the random and systematic components of the model (McCullagh and Nelder, 1989). Logistic regression belongs to the class of GLM. It involves fitting of a linear model to a function of the mean called the logit link function. Define y to be the response variable with y = 1 meaning a pixel was converted to residential land use between 1984 and 2000, and y = 0 otherwise. Then, y¼mþe
113 121,377 9 4
Standard deviation
(1)
where m is the mean of y and e is a vector of unobserved errors. The logit link g(m) is given by m gðmÞ ¼ ln (2) 1m and the logistic regression model is of the form m pðy ¼ 1jXÞ ¼ ln ¼ Xb gðmÞ ¼ ln 1m 1 pðy ¼ 1jXÞ (3) where p(y = 1jX) is the probability y takes the value 1 given the vector of independent variables X (Table 1), and b is a vector of unknown fixed effects with known model matrix X. The GLM assumes uncorrelated errors, so e = 0. Logistic regression was implemented using the
pˆ ðy ¼ 1jXÞ ¼
expbX ˆ 1 þ expbX
(4)
where the hat notation is used to indicate estimated values. 3.3. Combination of GLM with indicator kriging Indicator kriging (IK) is a parametric geostatistical tool for estimating the probability distribution of spatial variables on the basis of surrounding observations (Deutsch and Journel, 1998). It allows the use of binary variables. It may be used in estimating the probability of occurrence of an outcome into discrete classes. We used IK to estimate Pr[z(s0 = 1jz(s1), . . ., z(sn)], the conditional probability that a location s0 was converted to residential land use given the observations z(s1), . . ., z(sn), their locations, and their mutual dependence structure. The land use change at each location is a categorical variable coded in two mutually exclusive classes, ci c1 = 1, meaning a pixel is converted to residential land use between 1984 and 2000, and c2 = 0, otherwise. Indicator kriging took place in two steps. The first step was to model the spatial structure of each indicator variable using the semi-variogram. Secondly, the parameter estimates of the semi-variogram were used to predict the probability that a given location belongs to a given land use change class. Indicator kriging merely models the probability of change given the surrounding locations. A major interest, however, is to account for both the effects of the independent variables and the surrounding locations, on the probability of residential development. Thus, we combined logistic regression and indicator
442
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
kriging using Bayes’ theorem: pˆ ðy ¼ 1j½X [ sÞ ¼
pˆ 1 f ðXjc1 Þ pˆ 1 f ðXjc1 Þ þ pˆ 2 f ðXjc2 Þ
(5)
where pˆ ðy ¼ 1j½X [ sÞ is the probability that a location was converted to residential development given both the vector of independent variables X and the land use at location s, [ the algebraic union notation, pˆ 1 the estimated probability for residential land change, using indicator kriging, pˆ 2 the estimated probability ‘‘not converted to residential land use’’ using indicator kriging, and f(Xjci) is the conditional density of independent variables X given land use change class ci. Eq. (5) amounts to setting the estimated probabilities derived from IK as prior probabilities. As classification algorithms do not provide densities, we replaced f(Xjci) with the estimated probability using logistic regression. Goovaerts (2002) observed that when prior probabilities are equal for all classes as implemented in the GLM above, densities and posterior probabilities obtained from maximum likelihood estimates are proportional. Thus, replacing densities with probabilities estimated from logistic regression yields the same result as the theoretical expression in Eq. (5). 3.4. Generalized linear mixed model Generalized linear mixed models (GLMM) allow data from a variety of continuous and discrete distributions to be linked to a linear structure consisting of fixed and random effects. Like generalized linear models, GLMM assume normal random effects, and data can have discrete or continuous distribution in the exponential family (Breslow and Clayton, 1993). Correlation among the response variables are implicitly modeled in the GLMM procedure. The GLMM model is of the form gðmÞ ¼ Xb þ Zg þ e
(6)
where g is a vector of unknown random effects with known model matrix Z. The distribution of g is normal with E(g) = 0, and covariance matrix G. The distribution of e is also assumed to be normal with a covariance matrix R. Two approaches are normally employed in the GLMM procedure. The first approach models G by including random effects in the model specification. The second approach models correlations among the data directly by specifying a covariance structure for R. The second approach was adopted in this study to enable a direct comparison of estimated fixed effects with the GLM. The spatial
covariance was modeled as a parametric function of the distance dij between locations si and sj. The parameters of R are given by Ri j ¼ s 2 rðf; di j Þ
(7)
where s2 is the variance of the spatial process and r(f; dij) a correlation function with scale parameter f (Schabenberger and Pierce, 2002). The restricted maximum likelihood technique (REML) was used to simultaneously evaluate the fixed effects and covariance parameters. 3.5. Model comparison and validation Model validation was carried out in two fold. The first method uses the relative operating characteristic (ROC). The ROC compares binary data over the whole range of predicted probabilities. It aggregates into a single index of agreement, the ability of the model to predict the probability of residential land change at various locations on the landscape. That is, the ROC is a measure of the ability of the model to correctly specify location. An ROC curve plots the percentage of true-positive values versus the percentage of falsepositives in the model, where positivity is defined as a point being modeled as converted to residential land use. A true-positive value identifies a point that is classified as converted in both reality and by the model, whereas a false-positive value identifies a point that is classified as non-converted in reality and converted by the model. The area under the curve is the ROC statistic. It summarizes the ability of the model to discriminate between land use change classes. If the model assigns the probability of residential land change at random across the landscape, then the ROC will be equal to 0.5, whereas the ROC increases as the model assigns higher probabilities to points that are converted than points that are not converted. The second validation method uses the Wilcoxon non-parametric test to investigate differences in estimated probabilities for the models. The Wilcoxon test compares two or more groups of data. It makes no assumptions on the distributions of the data. It uses the ranks of the data rather than the raw values to calculate the test statistic. The null hypothesis is that the median is the same for each of the groups, whereas the alternative hypothesis is that the median is different in one of the groups. The test statistic was compared with critical values, which are based on group sample sizes.
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
4. Results and discussion 4.1. Spatial dependency of residential land change Semi-variograms describe the spatial dependency of spatial processes. Table 2 reveals substantial correlation in the independent variables. Four of the variables (population potential, number of different classes, distance from waterbodies and distance from waterworks) were fitted to exponential and Gaussian models that reach the sill asymptotically. Their practical range
443
varies from 6.1 km for number of different classes to 16.8 km for distance from waterworks. Frequency of non-urban land and travel time to the Central Business District were fitted to linear and power models, respectively, signifying continuous variability over the extent of the study area. Elevation and travel time to major roads have a sill, occurring at a range of 21.5 and 28.9 km, respectively. The fitted model for the indicator data is an exponential model (Table 2). The practical range of the indicator variogram is 16.8 km. This indicates that
Table 2 Parameters of the semi-variograms Nugget (C0) Independent variables ELEV 35 WATER 1.5 10 5 WATERWKS 2.0 10 5 CBD 0.29 ROADS 1.23 POP84 2.5 10 6 NURB 2.01 NDC 0.17 Dependent variables IK GLM GLMM
0.09 0.90 0.25
Partial sill (C1)
(Practical) range (km)
Slope
Exponent
Model
525 4.05 10 7 6.76 10 7 – 31.23 4.6 107 – 0.06
21.5 13 16.8 – 28.9 7.6 – 6.1
– – – 2.8 105 – – 4.50 –
– – – 1.67 – – – –
Spherical Gaussian Gaussian Power Spherical Exponential Linear Exponential
16.8 – 6.2
– 2.4 105 –
– – –
Exponential Linear Exponential
0.21 – 0.51
Table 3 Parameter estimates for the GLM and GLMM Parameter
Estimate
S.E.
Probability
95% confidence interval of estimates Lower limit
Upper limit
(a) GLM Intercept ELEV WATER NURB NDC POP84 CBD ROADS WATERWKS
8.545700 0.037200 0.000070 0.490000 0.782400 0.000032 0.031200 0.040100 0.000070
1.205500 0.003600 0.000014 0.101800 0.243500 8.95 106 0.004430 0.017400 0.000013
<0.0001 <0.0001 <0.0001 <0.0001 0.0013 0.0004 <0.0001 0.0213 <0.0001
– 0.030100 0.000100 0.689500 0.305100 0.000014 0.060000 0.074300 0.00009
– 0.04420 0.00005 0.29060 1.25970 0.000049 0.0427 0.00596 0.00004
(b) GLMM Intercept ELEV WATER NURB NDC POP84 CBD ROADS WATERWKS
5.262000 0.015700 0.000080 0.204000 0.195300 0.000022 0.031200 0.045740 0.000060
0.896100 0.003895 0.000027 0.046440 0.137500 0.000015 0.008768 0.018320 0.000029
<0.0001 <0.0001 0.0025 <0.0001 0.1557 0.1442 0.0004 0.0126 0.0514
– 0.008059 0.000140 0.295000 0.07437 7.55 106 0.048390 0.081670 0.000120
– 0.023340 0.000030 0.112900 0.465000 0.000052 0.014000 0.009810 3.50 107
444
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
there is correlation between change and non-change classes up to a separation distance of about 17 km. The spatially dependent component of variation constitutes about 72% of total variation. Such a long range confirms the existence of a strong trend within the study area. The micro-spatial variation (nugget) comprising 28% of the variation is indicative of a patchy terrain caused by the proliferation of illegal settlements on the landscape. The average low income earner has been priced out of formal land and housing market, leading to overcrowding and squatter settlements in several parts of Lagos. The spatial dependency of the GLM residuals was also modeled with a linear semi-variogram (Table 2), indicating continuous variability over the extent of the study area. This appears to be a reaction to the spatial variation of the independent variables. The nugget of the GLM (Pearson) residuals is 56% of the total variance. The total variance of 1.6 is higher than the theoretical value of 1, indicating that spatial correlation among sites still persists. Under this condition, the parameters of the GLM estimates are biased. The estimate of r, the spatial covariance parameter in R in the GLMM shows that spatial correlation between samples is negligible beyond a separation distance of 6.2 km (Table 2). This shorter range compared with that
of the indicator semi-variogram shows the ability of GLMM to account for interaction among sites. The estimated range corresponds to the average spatial extent of localities in Lagos. The total sill (C0 + C1) is lower than the theoretical value of 1, signifying the effectiveness of the GLMM in accounting for spatial correlation due to interaction among sites and reaction to trend. 4.2. Fixed effects for GLM and GLMM Table 3 shows that the signs of the fixed effects parameters are the same for the GLM and GLMM. Higher elevation favored residential land development due to higher costs associated with land reclamation in low lying areas. Residential land development more likely developed in areas with higher population potential in 1984. Conversion to residential land uses increased with increasing proximity to the Central Business District, major roads, publicly supplied water stations and waterbodies. The higher the proportion of non-urban land (forests and farmlands) in a neighborhood, the less likely the residential development. This is most likely due to the absence of infrastructure (electricity and pipe-water) in these areas. The more diversified the neighborhood, the higher the probability
Fig. 3. Actual conversion to residential land change and probability maps from models.
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
of residential land use conversion. All the independent variables are significant (P < 0.05) in the GLM, whereas NDC and POP84 are not significant in the GLMM. 4.3. Probability maps and validation of models
445
Table 4 ROC statistics to validate the models Model
Whole study area
Validation data (n = 588)
GLM GLM + IK GLMM
0.89 0.93 0.89
0.80 0.88 0.80
Maps depicting the probability of residential land development and the actual conversion between 1984 and 2000 are shown in Fig. 3. The three models generally assigned higher probabilities to locations that were truly converted than those that did not experience a change. Close inspection reveals differences in the spatial patterns of estimated probabilities. The GLM map shows noisy features (isolated pixels) towards the northeast and southeast of the study area. This stems from the lack of consideration for the spatial coordinates of the pixels in estimating probabilities. The GLM computes probability of change only in the parameter space of the independent variables. The GLMM and the GLM + IK on the other hand produced smoother maps. The GLM + IK further assigns very low probabilities to the northeast portion of the landscape that truly did not experience a change. The spatial agreement between modeled and actual residential land change as measured by the ROC is equal for GLM and GLMM, but higher for GLM + IK. Using smaller validation data (n = 588) slightly lowered the ROC statistic for the three models (Tables 4 and 5). The major benefits of the GLMM include providing quantitative information on the spatial dependency of land use change and reliability of its estimates. This is further demonstrated in Fig. 4, comparing the estimated probabilities of residential land change from marginal effects of the frequency of non-urban land for the GLM and GLMM. The GLM predicts lower and a more gradual increase in the probability of conversion to residential land use given a decrease
in the proportion of non-urban land in a neighborhood. Our observation during recent fieldwork in Lagos however shows that the GLMM estimates may be closer to reality. Given the escalating price of housing and urban land in the metropolis, there is aggressive extensification of residential areas through the conversion of cheaper non-urban land at the periphery of the metropolis. Application of estimates from the GLM could therefore lead to the wrong policy conclusions. The cumulative density functions of estimated probabilities for the validation data is shown in Fig. 5. The shape of the cdf for GLM + IK markedly differs from that of the other models. Its 25th percentile is 0.03 compared to 0.19 for GLMM and 0.20 for GLM. The estimated median probabilities are 0.72, 0.46 and 0.55 for GLM + IK, GLMM and GLM, respectively. The 75th percentile is located at 0.99 for GLM + IK, 0.67 for GLMM and 0.86 for GLM. Thus, the GLM + IK spreads predicted probabilities more evenly compared to the other models. This stems form the sensitivity of IK to nearby observations.
Fig. 4. Estimated change in probabilities of residential land development at different frequencies of non-urban land classes in 3 3 neighborhoods.
Fig. 5. Cumulative density functions for probability of residential land change.
Table 5 Wilcoxon test for estimated probabilities Models pair
Test statistic
Probability
GLM + IK vs. GLM GLMM vs. GLM GLMM vs. GLM + IK
2.46 12.52 8.32
0.014 <0.001 <0.001
446
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
5. Implications of the study and conclusions Spatial correlation contains useful information on the patterns of land use change. Variogram analysis as carried out in this study provides useful information for planners on the spatial determinants of land use. The continuous increase in the semi-variogram of travel time to the Central Business District may be due to spatial inequity in access to the major employment location. Spatial variation in travel time to major roads up to a distance of 29 km also suggests that localities are unequally served with major roads. This suggests the lack of integration of suburban areas to the metropolis and may hamper economic activities. This study demonstrates the benefits of geostatistical methods for improved understanding of the determinants of land use change. Variogram analysis proved useful in elucidating spatial correlation in land use data caused by trend and interaction among sites. A more accurate assessment of the impact of driving factors on the probability of land change was achieved by simultaneous analysis of correlation and fixed effects in the GLMM. The combination of GLM + IK may, in principle, be applied to large areas. Possibly, a Bayesian approach could be used, in which prior variogram estimation is done for data from one specific study area, and a Bayesian updating takes place for the new area (Cui et al., 1995). This leads to a Bayesian form of kriging, but avoids abundant data collection to have to fully estimate variogram parameters in a new region again. Acknowledgements Funding for this research was provided by the Japan Society for the Promotion of Science. The help
rendered by United Nations University Institute of Advanced Studies is also appreciated. We thank two reviewers for their comments on an earlier draft of the paper. References Breslow, N.E., Clayton, D.G., 1993. Approximate inference in generalized linear mixed models. J. Am. Stat. Assoc. 88, 9–25. Cui, H., Stein, A., Myers, D.E., 1995. Extension of spatial information, Bayesian kriging and updating of prior variogram parameters. Environmetrics 6, 373–384. Deutsch, C.V., Journel, A.G., 1998. Gslib Geostatistical Software Library and Users Guide. Oxford University Press, p. 369. Goovaerts, P., 2002. Geostatistical incorporation of spatial coordinates into supervised classification of hyperspectral data. J. Geograph. Syst. 4, 99–111. Griffith, D.A., Layne, L.J., 1999. A Casebook for Spatial Statistical Data Analysis. Oxford University Press, New York. Kruse, F.A., Boardman, J.W., Lefkoff, A.B., Heidebrecht, K.B., Shapiro, A.T., Barloon, P.J., Goetz, A.F.H., 1993. The spectral image processing system (SIPS)—interactive visualization and analysis of imaging spectrometer data. Remote Sens. Environ. 44, 145–163. Lambin, E.F., 1994. Modelling Deforestation Processes: A Review, TREES Publications Series B: Research Report No. 1, European Commission, EUR 15744 EN, 128 pp. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models, second ed. Chapman and Hall, New York. Overmars, K.P., de Koning, G.H.J., Veldkamp, A., 2003. Spatial autocorrelation in multi-scale land use change models. Ecol. Model. 164, 257–270. Schabenberger, O., Pierce, F.J., 2002. Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press, New York, p. 738. Vitousek, P.M., Mooney, H.A., Lubchenco, J., Melillo, J.M., 1997. Human domination of Earth’s ecosystems. Science 277, 494– 499. United Nations, 2002. World Urbanization Prospects: The 2001 Revision. United Nations, New York.
Geostatistical techniques for incorporating spatial correlation into land use change models Ademola K. Braimoh a,b,*, Takashi Onishi b a
United Nations University Institute of Advanced Studies, International Organizations Center, Yokohama 220-8502, Japan b Research Center for Advanced Science and Technology, University of Tokyo, 4-6-1 Komaba, 153-8904, Japan Received 21 March 2006; accepted 19 February 2007
Abstract Land use modeling requires large amounts of data that are typically spatially correlated. This study applies two geostatistical techniques to account for spatial correlation in residential land use change modeling. In the first approach, we combined generalized linear model (GLM) with indicator kriging to estimate the posterior probability of residential development. In the second approach, generalized linear mixed model (GLMM) was used to simultaneously model spatial correlation and regression fixed effects. Spatial agreement between actual and modeled land use change was higher for the GLM incorporating indicator kriging. The GLMM produced more reliable estimates and could be more useful in analyzing the effects of driving factors of land use change for land use planning. # 2007 Elsevier B.V. All rights reserved. Keywords: Generalized linear mixed model; Indicator kriging; Land use change; Lagos; Land use planning
1. Introduction Land use/cover change (LUCC) is one of the most important human impacts affecting ecosystems and environmental processes at local, regional and global scales. Such effects include land degradation and fragmentation, biodiversity loss, alteration in atmospheric composition and climate change (Vitousek et al., 1997). Major changes may yet be expected in the years to come as a result of population growth, and socio-economic and political developments. An understanding of the determinants of LUCC is therefore important for proper land use planning and the
* Corresponding author. Present address: Global Land Project, Sapporo Nodal Office, Hokkaido University, N9W8, Sapporo 060-0809, Japan. Tel.: +81 11 706 4531; fax: +81 11 706 4534. E-mail address: [email protected] (A.K. Braimoh). 0303-2434/$ – see front matter # 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jag.2007.02.005
development of sustainable environmental management policies. Models are required to formalize knowledge about land use change. They are important tools to explore the interactions between land use dynamics and the driving factors. By linking data and theory using formal equations, they help to improve our understanding of the processes underlying land use change. A land use model typically consists of three components: multitemporal land use maps derived from remotely sensed data, a multivariate function of the hypothesized determinants of change, and the resulting prediction map of land use change (Lambin, 1994). Land use models serve three main purposes. First, they help to identify the set of variables that are most important in explaining land use change. Second, they reveal the relative importance of the significant variables. The magnitude of the effect of factors on land change is important because it may indicate areas where policy
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
intervention is most urgently needed, and whether local action to mitigate adverse effect of land change will be sustainable. Last, land use models may be useful in predicting the location of land use change in the immediate future. The ability to predict vulnerable places to land use change helps to mitigate the adverse effects of land change. Land use modeling often requires biophysical and social data that are typically spatially correlated. Spatial correlation is the tendency of objects near in space to be more similar than those further apart. It is a nuisance in conventional statistics, but in land use modeling, it
439
reveals the pattern, and possibly the processes associated with the pattern of land use. Spatial correlation in land use derives from two sources (Overmars et al., 2003). The first source is the distribution of the biophysical and social variables that determines land use pattern and change. The second source is spatial proximity (interaction) between land use types. Spatial interaction in land use might be due to individual choices in land use location, or due to spatial policies in the allocation of land. For instance, the proximity to publicly supplied infrastructure, employment locations or the need for social interaction might
Fig. 1. The study area.
440
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
encourage the clustering of residential land use in an area. Similarly, spatial planning policies may encourage the concentration of business premises in designated commercial centers. Conversely, down-zoning policies may restrict the extent and intensity of urban development in an area. It is therefore essential to account for spatial correlation in land use change models. Ignoring spatial correlation in land use modeling may result in false detection of significant relationships and biased model parameter estimates (Griffith and Layne, 1999). The objectives of this study are to apply and compare two geostatistical techniques for incorporating modeling residential land change. The first method combines logistic regression with Indicator kriging (IK), whereas the second method concurrently models spatial correlation and regression fixed effects. The study applies to residential land use change in Lagos, Nigeria. 2. The study area The study area of approximately 2910 km2 lies between latitudes 68260 and 68500 N and stretches between longitudes 38090 and 38460 E (Fig. 1). Lagos is comprised of depositional landform features: wetlands, barrier islands, beaches, low-lying tidal flats and estuaries. It is characterized by a wet equatorial climate with mean annual rainfall above 1800 mm, mean monthly temperature of 30 8C, and relative humidity ranging from 80 to 100%. The terrain is relatively flat with a mean elevation of about 24 m. The physical nature of the terrain determines the pattern of settlement, and markedly limits the supply of land for urban development. Lagos is the dominant center for non-agricultural production, distribution and business services in Nigeria. Its current population of over 13 million is estimated to grow annually at the rate of 5.6% (UN, 2002). Housing is grossly inadequate to cater for the increasing population. This problem is related to a number of social and institutional factors such as the skewed distribution of private land, the high costs of undeveloped land and housing construction, and weak planning regulations. Thus, new settlements, even if substandard, are only affordable and are springing up at the periphery of the Lagos metropolis. Despite the importance of Lagos in the socio-economic development of Nigeria and for the West African sub-region, there is limited study of the factors responsible for the development of its spatial pattern of land use. Such a dearth of knowledge hinders the development of policies for sustainable land use planning and management of the city.
3. Methods 3.1. Data sources Land use/cover maps of the study area for 1984 and 2000 prepared from Landsat TM images were used. Land cover classification was performed using spectral angle mapping technique (Kruse et al., 1993). Spectral endmembers for the classification were derived from existing land use maps and intensive field survey in May 2005. Fig. 2 shows the maps for four land classes: residential, industrial/commercial, non-urban land, and water; and their relative proportions in 1984 and 2000. The dependent variable for this study, residential development between 1984 and 2000, was extracted from the maps within a geographical information system. Eight spatially explicit independent variables hypothesized to affect residential development were also developed within a GIS (Table 1). These include travel times to major roads and the central business district, Euclidean distances from waterbodies and publicly supplied water stations, population potential in 1984, elevation, frequency of non-urban land within 3 3 neighborhoods in 1984, and number of different land classes within 3 3 neighborhoods in 1984. All variables were mapped at a resolution of 90 m, dictated by the grain size of the digital elevation model. A total
Fig. 2. Land use/cover maps and proportions.
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
441
Table 1 Independent variables list Variable name/description
Abbreviation
Minimum
Mean
Maximum
Elevation (m) Distance from rivers and other waterbodies (m) Distance from waterworks (m) Travel time to Lagos Island, the Central Business District (min) Travel time to major roads (min) Population potential in 1984 Number of non-urban land within 270 m neighborhood Number of different classes within 270 m neighborhood
ELEV WATER WATERWKS CBD
36 0 0 0
24 6,665 13,026 107
115 26,364 36,495 179
ROADS POP84 NURB NDC
0 334 0 1
of 2668 points were randomly selected, out of which 2080 points were used for modeling and the remaining 588 points for validation.
56 6,404 5.99 1.27
23 6532 8337 32 33 8046 3.94 0.47
Coefficient of variation (%) 96 98 64 30 59 126 66 37
maximum likelihood algorithm with equal prior probabilities. After back transformation, the result of the logistic regression may be expressed in terms of conditional probability as
3.2. Generalized linear model (GLM)
ˆ
Generalized Linear Model relaxes the independent and identically distributed assumption commonly used for normal data. It has three distinct parts: a random component which is the qualitative variable being modeled, a systematic or fixed effect component which refers to the independent variables that could be either categorical or continuous; and a link function which is used to relate the random and systematic components of the model (McCullagh and Nelder, 1989). Logistic regression belongs to the class of GLM. It involves fitting of a linear model to a function of the mean called the logit link function. Define y to be the response variable with y = 1 meaning a pixel was converted to residential land use between 1984 and 2000, and y = 0 otherwise. Then, y¼mþe
113 121,377 9 4
Standard deviation
(1)
where m is the mean of y and e is a vector of unobserved errors. The logit link g(m) is given by m gðmÞ ¼ ln (2) 1m and the logistic regression model is of the form m pðy ¼ 1jXÞ ¼ ln ¼ Xb gðmÞ ¼ ln 1m 1 pðy ¼ 1jXÞ (3) where p(y = 1jX) is the probability y takes the value 1 given the vector of independent variables X (Table 1), and b is a vector of unknown fixed effects with known model matrix X. The GLM assumes uncorrelated errors, so e = 0. Logistic regression was implemented using the
pˆ ðy ¼ 1jXÞ ¼
expbX ˆ 1 þ expbX
(4)
where the hat notation is used to indicate estimated values. 3.3. Combination of GLM with indicator kriging Indicator kriging (IK) is a parametric geostatistical tool for estimating the probability distribution of spatial variables on the basis of surrounding observations (Deutsch and Journel, 1998). It allows the use of binary variables. It may be used in estimating the probability of occurrence of an outcome into discrete classes. We used IK to estimate Pr[z(s0 = 1jz(s1), . . ., z(sn)], the conditional probability that a location s0 was converted to residential land use given the observations z(s1), . . ., z(sn), their locations, and their mutual dependence structure. The land use change at each location is a categorical variable coded in two mutually exclusive classes, ci c1 = 1, meaning a pixel is converted to residential land use between 1984 and 2000, and c2 = 0, otherwise. Indicator kriging took place in two steps. The first step was to model the spatial structure of each indicator variable using the semi-variogram. Secondly, the parameter estimates of the semi-variogram were used to predict the probability that a given location belongs to a given land use change class. Indicator kriging merely models the probability of change given the surrounding locations. A major interest, however, is to account for both the effects of the independent variables and the surrounding locations, on the probability of residential development. Thus, we combined logistic regression and indicator
442
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
kriging using Bayes’ theorem: pˆ ðy ¼ 1j½X [ sÞ ¼
pˆ 1 f ðXjc1 Þ pˆ 1 f ðXjc1 Þ þ pˆ 2 f ðXjc2 Þ
(5)
where pˆ ðy ¼ 1j½X [ sÞ is the probability that a location was converted to residential development given both the vector of independent variables X and the land use at location s, [ the algebraic union notation, pˆ 1 the estimated probability for residential land change, using indicator kriging, pˆ 2 the estimated probability ‘‘not converted to residential land use’’ using indicator kriging, and f(Xjci) is the conditional density of independent variables X given land use change class ci. Eq. (5) amounts to setting the estimated probabilities derived from IK as prior probabilities. As classification algorithms do not provide densities, we replaced f(Xjci) with the estimated probability using logistic regression. Goovaerts (2002) observed that when prior probabilities are equal for all classes as implemented in the GLM above, densities and posterior probabilities obtained from maximum likelihood estimates are proportional. Thus, replacing densities with probabilities estimated from logistic regression yields the same result as the theoretical expression in Eq. (5). 3.4. Generalized linear mixed model Generalized linear mixed models (GLMM) allow data from a variety of continuous and discrete distributions to be linked to a linear structure consisting of fixed and random effects. Like generalized linear models, GLMM assume normal random effects, and data can have discrete or continuous distribution in the exponential family (Breslow and Clayton, 1993). Correlation among the response variables are implicitly modeled in the GLMM procedure. The GLMM model is of the form gðmÞ ¼ Xb þ Zg þ e
(6)
where g is a vector of unknown random effects with known model matrix Z. The distribution of g is normal with E(g) = 0, and covariance matrix G. The distribution of e is also assumed to be normal with a covariance matrix R. Two approaches are normally employed in the GLMM procedure. The first approach models G by including random effects in the model specification. The second approach models correlations among the data directly by specifying a covariance structure for R. The second approach was adopted in this study to enable a direct comparison of estimated fixed effects with the GLM. The spatial
covariance was modeled as a parametric function of the distance dij between locations si and sj. The parameters of R are given by Ri j ¼ s 2 rðf; di j Þ
(7)
where s2 is the variance of the spatial process and r(f; dij) a correlation function with scale parameter f (Schabenberger and Pierce, 2002). The restricted maximum likelihood technique (REML) was used to simultaneously evaluate the fixed effects and covariance parameters. 3.5. Model comparison and validation Model validation was carried out in two fold. The first method uses the relative operating characteristic (ROC). The ROC compares binary data over the whole range of predicted probabilities. It aggregates into a single index of agreement, the ability of the model to predict the probability of residential land change at various locations on the landscape. That is, the ROC is a measure of the ability of the model to correctly specify location. An ROC curve plots the percentage of true-positive values versus the percentage of falsepositives in the model, where positivity is defined as a point being modeled as converted to residential land use. A true-positive value identifies a point that is classified as converted in both reality and by the model, whereas a false-positive value identifies a point that is classified as non-converted in reality and converted by the model. The area under the curve is the ROC statistic. It summarizes the ability of the model to discriminate between land use change classes. If the model assigns the probability of residential land change at random across the landscape, then the ROC will be equal to 0.5, whereas the ROC increases as the model assigns higher probabilities to points that are converted than points that are not converted. The second validation method uses the Wilcoxon non-parametric test to investigate differences in estimated probabilities for the models. The Wilcoxon test compares two or more groups of data. It makes no assumptions on the distributions of the data. It uses the ranks of the data rather than the raw values to calculate the test statistic. The null hypothesis is that the median is the same for each of the groups, whereas the alternative hypothesis is that the median is different in one of the groups. The test statistic was compared with critical values, which are based on group sample sizes.
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
4. Results and discussion 4.1. Spatial dependency of residential land change Semi-variograms describe the spatial dependency of spatial processes. Table 2 reveals substantial correlation in the independent variables. Four of the variables (population potential, number of different classes, distance from waterbodies and distance from waterworks) were fitted to exponential and Gaussian models that reach the sill asymptotically. Their practical range
443
varies from 6.1 km for number of different classes to 16.8 km for distance from waterworks. Frequency of non-urban land and travel time to the Central Business District were fitted to linear and power models, respectively, signifying continuous variability over the extent of the study area. Elevation and travel time to major roads have a sill, occurring at a range of 21.5 and 28.9 km, respectively. The fitted model for the indicator data is an exponential model (Table 2). The practical range of the indicator variogram is 16.8 km. This indicates that
Table 2 Parameters of the semi-variograms Nugget (C0) Independent variables ELEV 35 WATER 1.5 10 5 WATERWKS 2.0 10 5 CBD 0.29 ROADS 1.23 POP84 2.5 10 6 NURB 2.01 NDC 0.17 Dependent variables IK GLM GLMM
0.09 0.90 0.25
Partial sill (C1)
(Practical) range (km)
Slope
Exponent
Model
525 4.05 10 7 6.76 10 7 – 31.23 4.6 107 – 0.06
21.5 13 16.8 – 28.9 7.6 – 6.1
– – – 2.8 105 – – 4.50 –
– – – 1.67 – – – –
Spherical Gaussian Gaussian Power Spherical Exponential Linear Exponential
16.8 – 6.2
– 2.4 105 –
– – –
Exponential Linear Exponential
0.21 – 0.51
Table 3 Parameter estimates for the GLM and GLMM Parameter
Estimate
S.E.
Probability
95% confidence interval of estimates Lower limit
Upper limit
(a) GLM Intercept ELEV WATER NURB NDC POP84 CBD ROADS WATERWKS
8.545700 0.037200 0.000070 0.490000 0.782400 0.000032 0.031200 0.040100 0.000070
1.205500 0.003600 0.000014 0.101800 0.243500 8.95 106 0.004430 0.017400 0.000013
<0.0001 <0.0001 <0.0001 <0.0001 0.0013 0.0004 <0.0001 0.0213 <0.0001
– 0.030100 0.000100 0.689500 0.305100 0.000014 0.060000 0.074300 0.00009
– 0.04420 0.00005 0.29060 1.25970 0.000049 0.0427 0.00596 0.00004
(b) GLMM Intercept ELEV WATER NURB NDC POP84 CBD ROADS WATERWKS
5.262000 0.015700 0.000080 0.204000 0.195300 0.000022 0.031200 0.045740 0.000060
0.896100 0.003895 0.000027 0.046440 0.137500 0.000015 0.008768 0.018320 0.000029
<0.0001 <0.0001 0.0025 <0.0001 0.1557 0.1442 0.0004 0.0126 0.0514
– 0.008059 0.000140 0.295000 0.07437 7.55 106 0.048390 0.081670 0.000120
– 0.023340 0.000030 0.112900 0.465000 0.000052 0.014000 0.009810 3.50 107
444
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
there is correlation between change and non-change classes up to a separation distance of about 17 km. The spatially dependent component of variation constitutes about 72% of total variation. Such a long range confirms the existence of a strong trend within the study area. The micro-spatial variation (nugget) comprising 28% of the variation is indicative of a patchy terrain caused by the proliferation of illegal settlements on the landscape. The average low income earner has been priced out of formal land and housing market, leading to overcrowding and squatter settlements in several parts of Lagos. The spatial dependency of the GLM residuals was also modeled with a linear semi-variogram (Table 2), indicating continuous variability over the extent of the study area. This appears to be a reaction to the spatial variation of the independent variables. The nugget of the GLM (Pearson) residuals is 56% of the total variance. The total variance of 1.6 is higher than the theoretical value of 1, indicating that spatial correlation among sites still persists. Under this condition, the parameters of the GLM estimates are biased. The estimate of r, the spatial covariance parameter in R in the GLMM shows that spatial correlation between samples is negligible beyond a separation distance of 6.2 km (Table 2). This shorter range compared with that
of the indicator semi-variogram shows the ability of GLMM to account for interaction among sites. The estimated range corresponds to the average spatial extent of localities in Lagos. The total sill (C0 + C1) is lower than the theoretical value of 1, signifying the effectiveness of the GLMM in accounting for spatial correlation due to interaction among sites and reaction to trend. 4.2. Fixed effects for GLM and GLMM Table 3 shows that the signs of the fixed effects parameters are the same for the GLM and GLMM. Higher elevation favored residential land development due to higher costs associated with land reclamation in low lying areas. Residential land development more likely developed in areas with higher population potential in 1984. Conversion to residential land uses increased with increasing proximity to the Central Business District, major roads, publicly supplied water stations and waterbodies. The higher the proportion of non-urban land (forests and farmlands) in a neighborhood, the less likely the residential development. This is most likely due to the absence of infrastructure (electricity and pipe-water) in these areas. The more diversified the neighborhood, the higher the probability
Fig. 3. Actual conversion to residential land change and probability maps from models.
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
of residential land use conversion. All the independent variables are significant (P < 0.05) in the GLM, whereas NDC and POP84 are not significant in the GLMM. 4.3. Probability maps and validation of models
445
Table 4 ROC statistics to validate the models Model
Whole study area
Validation data (n = 588)
GLM GLM + IK GLMM
0.89 0.93 0.89
0.80 0.88 0.80
Maps depicting the probability of residential land development and the actual conversion between 1984 and 2000 are shown in Fig. 3. The three models generally assigned higher probabilities to locations that were truly converted than those that did not experience a change. Close inspection reveals differences in the spatial patterns of estimated probabilities. The GLM map shows noisy features (isolated pixels) towards the northeast and southeast of the study area. This stems from the lack of consideration for the spatial coordinates of the pixels in estimating probabilities. The GLM computes probability of change only in the parameter space of the independent variables. The GLMM and the GLM + IK on the other hand produced smoother maps. The GLM + IK further assigns very low probabilities to the northeast portion of the landscape that truly did not experience a change. The spatial agreement between modeled and actual residential land change as measured by the ROC is equal for GLM and GLMM, but higher for GLM + IK. Using smaller validation data (n = 588) slightly lowered the ROC statistic for the three models (Tables 4 and 5). The major benefits of the GLMM include providing quantitative information on the spatial dependency of land use change and reliability of its estimates. This is further demonstrated in Fig. 4, comparing the estimated probabilities of residential land change from marginal effects of the frequency of non-urban land for the GLM and GLMM. The GLM predicts lower and a more gradual increase in the probability of conversion to residential land use given a decrease
in the proportion of non-urban land in a neighborhood. Our observation during recent fieldwork in Lagos however shows that the GLMM estimates may be closer to reality. Given the escalating price of housing and urban land in the metropolis, there is aggressive extensification of residential areas through the conversion of cheaper non-urban land at the periphery of the metropolis. Application of estimates from the GLM could therefore lead to the wrong policy conclusions. The cumulative density functions of estimated probabilities for the validation data is shown in Fig. 5. The shape of the cdf for GLM + IK markedly differs from that of the other models. Its 25th percentile is 0.03 compared to 0.19 for GLMM and 0.20 for GLM. The estimated median probabilities are 0.72, 0.46 and 0.55 for GLM + IK, GLMM and GLM, respectively. The 75th percentile is located at 0.99 for GLM + IK, 0.67 for GLMM and 0.86 for GLM. Thus, the GLM + IK spreads predicted probabilities more evenly compared to the other models. This stems form the sensitivity of IK to nearby observations.
Fig. 4. Estimated change in probabilities of residential land development at different frequencies of non-urban land classes in 3 3 neighborhoods.
Fig. 5. Cumulative density functions for probability of residential land change.
Table 5 Wilcoxon test for estimated probabilities Models pair
Test statistic
Probability
GLM + IK vs. GLM GLMM vs. GLM GLMM vs. GLM + IK
2.46 12.52 8.32
0.014 <0.001 <0.001
446
A.K. Braimoh, T. Onishi / International Journal of Applied Earth Observation and Geoinformation 9 (2007) 438–446
5. Implications of the study and conclusions Spatial correlation contains useful information on the patterns of land use change. Variogram analysis as carried out in this study provides useful information for planners on the spatial determinants of land use. The continuous increase in the semi-variogram of travel time to the Central Business District may be due to spatial inequity in access to the major employment location. Spatial variation in travel time to major roads up to a distance of 29 km also suggests that localities are unequally served with major roads. This suggests the lack of integration of suburban areas to the metropolis and may hamper economic activities. This study demonstrates the benefits of geostatistical methods for improved understanding of the determinants of land use change. Variogram analysis proved useful in elucidating spatial correlation in land use data caused by trend and interaction among sites. A more accurate assessment of the impact of driving factors on the probability of land change was achieved by simultaneous analysis of correlation and fixed effects in the GLMM. The combination of GLM + IK may, in principle, be applied to large areas. Possibly, a Bayesian approach could be used, in which prior variogram estimation is done for data from one specific study area, and a Bayesian updating takes place for the new area (Cui et al., 1995). This leads to a Bayesian form of kriging, but avoids abundant data collection to have to fully estimate variogram parameters in a new region again. Acknowledgements Funding for this research was provided by the Japan Society for the Promotion of Science. The help
rendered by United Nations University Institute of Advanced Studies is also appreciated. We thank two reviewers for their comments on an earlier draft of the paper. References Breslow, N.E., Clayton, D.G., 1993. Approximate inference in generalized linear mixed models. J. Am. Stat. Assoc. 88, 9–25. Cui, H., Stein, A., Myers, D.E., 1995. Extension of spatial information, Bayesian kriging and updating of prior variogram parameters. Environmetrics 6, 373–384. Deutsch, C.V., Journel, A.G., 1998. Gslib Geostatistical Software Library and Users Guide. Oxford University Press, p. 369. Goovaerts, P., 2002. Geostatistical incorporation of spatial coordinates into supervised classification of hyperspectral data. J. Geograph. Syst. 4, 99–111. Griffith, D.A., Layne, L.J., 1999. A Casebook for Spatial Statistical Data Analysis. Oxford University Press, New York. Kruse, F.A., Boardman, J.W., Lefkoff, A.B., Heidebrecht, K.B., Shapiro, A.T., Barloon, P.J., Goetz, A.F.H., 1993. The spectral image processing system (SIPS)—interactive visualization and analysis of imaging spectrometer data. Remote Sens. Environ. 44, 145–163. Lambin, E.F., 1994. Modelling Deforestation Processes: A Review, TREES Publications Series B: Research Report No. 1, European Commission, EUR 15744 EN, 128 pp. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models, second ed. Chapman and Hall, New York. Overmars, K.P., de Koning, G.H.J., Veldkamp, A., 2003. Spatial autocorrelation in multi-scale land use change models. Ecol. Model. 164, 257–270. Schabenberger, O., Pierce, F.J., 2002. Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press, New York, p. 738. Vitousek, P.M., Mooney, H.A., Lubchenco, J., Melillo, J.M., 1997. Human domination of Earth’s ecosystems. Science 277, 494– 499. United Nations, 2002. World Urbanization Prospects: The 2001 Revision. United Nations, New York.