A spatial multinomial logit model for analysing urban expansion
aa r X i v : . [ ec on . E M ] A ug A spatial multinomial logit model for analysing urbanexpansion
Tamás Krisztin ∗1 , Philipp Piribauer , and Michael Wögerer International Institute for Applied Systems Analysis (IIASA) Austrian Institute of Economic Research (WIFO)
Abstract
The paper proposes a Bayesian multinomial logit model to analyse spatial patterns of urbanexpansion. The specification assumes that the log-odds of each class follow a spatial au-toregressive process. Using recent advances in Bayesian computing, our model allows fora computationally efficient treatment of the spatial multinomial logit model. This allows usto assess spillovers between regions and across land use classes. In a series of Monte Carlostudies, we benchmark our model against other competing specifications. The paper alsoshowcases the performance of the proposed specification using European regional data. Ourresults indicate that spatial dependence plays a key role in land sealing process of croplandand grassland. Moreover, we uncover land sealing spillovers across multiple classes of arableland.
Keywords: urban expansion, land use change, spatial multinomial logit model, Europeanregions
JEL Codes:
C11, C21, C25, O13, R14 ∗ Corresponding author : Tamás Krisztin, International Institute for Applied Systems Analysis (IIASA), Schloßplatz1, 2361 Laxenburg, Austria.
E-mail : [email protected]. The research carried out in this paper was supported by fundsof the Oesterreichische Nationalbank (Jubilaeumsfond project number: 18043). Introduction
Increased urbanisation and expansion of cities as a direct result of economic and population growth,coupled with intensifying climate change, poses a key challenge for policy makers (IPBES, 2019).The location choice of new urban developments is of particular importance, as land is a finiteresource. Expanding artificial surfaces is both expensive and time consuming to reverse, resultingin long-term impacts on land use and land cover. The conversion of natural habitats to artificialsurfaces thus has a direct and potentially irreversible impact on biodiversity (Leclère et al., 2018).On the other hand, if arable land is built up, global food security is threatened and urban expansionmight spillover to other types of land use.Conversion of land to urban surfaces is a decision usually taken by the landowners, which areeither regional governments or private land-holders. In an economic framework, this decision isunderstood as a trade-off between the relative profitabilities of land uses and respective conversioncosts (Miller and Plantinga, 1999). Potential profits from land ownership are typically assessedusing various proxies for land rents (Chakir and Lungarska, 2017), while conversion costs relyon the quality of land and national regulations restricting land transformation. Land use changemodels targeting aggregate administrative levels focus on capturing the outcomes of regionalpolicies (Ay et al., 2017). This is of special importance within the EU, where regional policies(such as the structural funds) are aimed at this level (Alexiadis et al., 2013).Within a regional econometric framework, the land use expansion decision can be modeled asa random choice, with the multinomial logit model representing a popular option (Lubowski et al.,2008; Chakir, 2009). The particular advantage is that a joint modelling of land sealing processescan take into account spillovers across land use classes. When dealing with compositional (shares)data for land use, the multinomial logit random choice model can either be estimated directly fromthe multinomial logit form (Li et al., 2013) or from its log-linearized form (Chakir and Lungarska,2017). While the log-linearized version of the model represents a popular choice due to its easeof transformation, it suffers from the usual problems of log-transformation, namely that frequentlyland use shares are zero and accommodating these observations inherently biases the estimates.Spatial dependence, both from unobserved spatially varying variables, as well as contin-gent on the choice of neighbouring regions, is well documented in the land use choice litera-ture (Chakir and Parent, 2009; Chakir and Le Gallo, 2013; Li et al., 2013). In a regional econo-metrics context, a wide number of studies stress the inherent importance of spatial spillovers(LeSage and Pace, 2009). When estimating models for land use change on a small-scale level, theproblem of spatial dependence becomes even more central. Specifically, neglecting to account forspatial autocorrelation may result in severely biased estimates and erroneous policy conclusions.However, spatial dependence in multinomial logit frameworks has been so far neglected by thespatial econometric literature, with the exception of GMM based approaches (Klier and McMillen,2008).Within this paper our contribution to the existing literature is twofold. First and foremost wepresent a novel Bayesian approach for capturing spatial dependence among land use changes usinga multinomial logit framework. By combining the spatial autoregressive and multionomial logit2rameworks, our specification can account for cross-regional and cross-land use class spillovers.The estimation approach builds on recent advances in Bayesian modelling of logit type models(Krisztin and Piribauer, 2020) and employs latent Pólya-Gamma distributed variables. We demon-strate the virtues of our approach in a series of Monte Carlo studies. Our second contribution is anovel examination land use change processes on a regional pan-European level. For this we relyon an extensive dataset of land use changes to assess the share of urban gains originating fromcropland, grassland, forest and other fallow land. Our framework allows us to shed light on thesmall-scale spatial dynamics of land sealing processes in European regions.The remainder of the paper is structured as follows. Section 2 outlines the theoretical model ofurban expansion, as well as its multinomial logit variant. Section 3 focuses on the estimation frame-work. Section 4 presents the Monte Carlo benchmarks of the proposed econometric estimationapproach. Section 5 presents the results for urban expansion in Europe. Section 6 concludes.
In this paper, we estimate an econometric model which aims at explaining the choice of land buyers(both public and private) for the purpose of converting it to urban, artificial surfaces in N regions.In a given region i (with i = , ..., N ), land buyers may acquire land from J different land uses. Inour case these are cropland, grassland, forest, and other natural land. Within a region the buyersare assumed to be price taker and their choices are assumed to be homogeneous and risk-neutral.In an economic sense this constitutes a profit maximization problem of land buyers (Lubowski et al.,2008; Miller and Plantinga, 1999), which is directly dependent on the associated profits and costsof the converted land. In addition, to account for the expected net present value of rents fromurban land use and the respective conversion costs, land buyers also face the opportunity costs ofalternatives usages.Such frameworks have been adopted among others by Lubowski et al. (2008), Chakir and Parent(2009), and Li et al. (2013). For estimation of parameters relating to observed buyers’ choices, theprofit maximization problem can be formulated within a multinomial limited dependent variableframework. Let y ij be the observed share of urban expansion from land use j relative to thetotal urban expansion in region i . Econometric estimation thus concerns itself with modellingthe probability of observing y ij . Within the multinomial logit framework, this probability can bemodelled as a function of choice specific log-odds µ ij , weighted by the sum of log-odds over all When land use specific observations y ij are shares, a popular choice for estimating the multinomial model isto apply a log-linear transformation, where the dependent variables correspond to log ( y ij / y iJ ) and perform standardregression analysis (see, for example Chakir and Lungarska, 2017; Chakir, 2009). The main drawbacks of this approachare twofold. First, Jensen’s inequality states that the expectation of a logarithm is not equal to the logarithm of theexpectation. Therefore, log-linearization inherently introduces a bias in the estimated slope coefficients. Second, inempirical applications frequently a large number of observed choices y ij are equal to zero, thus necessitating either acensoring of observations or adding a constant to all observations, both of which have been demonstrated to lead tosubstantial bias. µ ij ′ ( j ′ = , . . . , J ): p ( y ij ) = exp µ ij Í Jj ′ = exp µ ij ′ . (2.1)In the standard non-spatial multinomial framework µ ij is specified as a function of k explanatoryvariables, with corresponding choice-specific slope coefficients, which are to be estimated. Theexplanatory variables correspond to the expected rents and conversion costs with respect to landuse j .Spatial dependence among log-odds µ ij in Eq. (2.1) involves the assumption that the choices ofurban land buyers do not solely depend on rent and conversion costs in their own region i , but alsoon other regions’ characteristics as well. This assumption implies that the probability of observinga land use choice in region i also depends on land use choices of all other regions. This assumptionis based on the spatial nature of land expansion: before construction, investors typically scopemultiple investment opportunities, which might not be contiguous, but located across regions inspatial proximity to each other.Following the spatial econometric literature, such dependencies can be incorporated by im-posing an exogenous neighbourhood structure through a non-negative and row-stochastic spatialweight matrix. Let W be such an N × N spatial weight matrix. Two regions i and i ′ are assumedto be neighbours of w ii ′ >
0, otherwise w ii ′ =
0. No region is a neighbour to itself, thus w ii = µ j = ρ j W µ j + X β j + ε j µ j = A − ( X β j + ε j ) , (2.2)with A − j = ( I N − ρ j W ) − where I N denotes an N × N identity matrix. The N × K matrix X = [ x , . . . , x K ] collects the K vectors of explanatory variables and β j denote the respective K × j . The N × ε j contains independentlyand identically Gaussian distributed disturbance terms, with zero mean and σ j variance. The(scalar) parameter ρ j measures the strength of spatial autocorrelation for land use class j , withsufficient stability condition ρ j ∈ (− , ) , where positive (negative) values of ρ indicate positive(negative) spatial autocorrelation. Note, that the model allows for different ρ j across land useclasses. In the absence of spatial autocorrelation ( ρ = ... = ρ J = N × µ j = [ µ j , . . . , µ N j ] ′ thus also depend on the characteristics of other regions in the sample.Spatial dependence is introduced by the spatial multiplier A − j = ( I N − ρ j W ) − = Í ∞ r = ρ rj W r . With the identifying restriction that the spatial autocorrelation coefficient associated with the J -th land use class ρ J = It is worth noting that the standard SAR model can be extended to more flexible spatial econometric modelspecifications in a straightforward way. Specifically, one may additionally include spatially lagged explanatory variables,resulting in a so-called spatial Durbin model (SDM) specification (see, for example LeSage and Pace 2009). A similarextension is presented in the empirical exercise. i , would not only result in changes of the observed shares y ij inthe own region, but in other regions as well. Through the nature of the multinomial logit model,where marginal impacts to one choice j also affect the shares of all other choices, this impliesthat in a spatial dependent setting marginal impacts of y ij have spillover effects over regions andchoices as well.Note, that through the normally distributed residual error vector ε j , the residuals of the modelin Eq. (2.2) are effectively decomposed into two components: first, the heteroscedastic errors aris-ing from the logistic model in Eq. (2.1) and second, the normally distributed error term ε j with σ j variance. Spatial dependence in the errors is captured through the latter. Similar to spatial autore-gressive variants of standard probit (LeSage et al., 2011) and logit models (Krisztin and Piribauer2020), innovation variances σ j are restricted to unity, in order to identify the logistic errors. We propose a Bayesian estimation strategy for the SAR multinomial logit model, which builds on theidea of introducing a latent variable in order to facilitate the estimation of the multinomial logit like-lihood. This estimation strategy has been widely employed in recent Bayesian econometric litera-ture for tackling models featuring non Gaussian distributions (see e.g. Frühwirth-Schnatter and Frühwirth,2012; Frühwirth-Schnatter et al., 2009). To illustrate the core problem, consider the likelihood ofthe multinomial logit model in Eq. (2.1): N Ö i = J Ö j = (cid:0) exp µ ij (cid:1) y i j Í Jj ′ = exp µ ij ′ . (3.1)Note, that the likelihood contribution of observation i relies not only on µ ij , but on the log-odds ofmaking other choices as well. This well-known non-linearity in the likelihood greatly complicatesthe estimation of the unknown slope and spatial autoregressive coefficients.Within a Bayesian framework the focus of estimation frequently lies mainly on finding condi-tional posterior distributions for the parameters of interest. In fact, assuming suitable priors p ( β j ) ,the conditional posterior of β j can be expressed conditional on all other slope coefficients β − j and ρ (see Holmes and Held, 2006): p (cid:16) β j | β − j , ρ (cid:17) = p (cid:16) β j (cid:17) N Ö i = (cid:18) exp η ij + exp η ij (cid:19) y i j (cid:18) + exp η ij (cid:19) − y i j (3.2)with η ij = µ ij − C ij and C ij = log J Õ j ′ , j exp µ ij ′ . While this distribution cannot be easily sampled from, we follow the work of Polson et al. (2013),which has been adopted to the spatial autoregressive variant of a bivariate logit distribution(Krisztin and Piribauer, 2020). A particularly useful result in Polson et al. (2013) is the fact thatconditional on introducing a Pólya-Gamma distributed latent random variable, exponential type5istributions such as the one in Eq. (3.2) can be recast as Gaussian, where posterior sampling canbe easily achieved.Particularly, when conditioning on ω ij ∼ PG( , ) – where PG( , ) denotes a Pólya-Gammadistribution with rate one and shape zero – the conditional posterior of the slope parametersassociated with choice j can be reformulated as: p (cid:16) β j | β − j , ρ, ω j (cid:17) ∝ p (cid:16) β j (cid:17) N Ö i = exp (cid:0) κ ij η ij (cid:1) exp η ij ω ij ! PG( ω ij | , )∝ p (cid:16) β j (cid:17) exp (cid:26) − (cid:16) (cid:2) z j − c j (cid:3) − A − j X (cid:17) ′ Ω j (cid:16) (cid:2) z j − c j (cid:3) − A − j X (cid:17) (cid:27) , where ω j = [ ω j , ..., ω N j ] ′ and κ ij = y ij − /
2. The conditional posterior has working responses z j = [ κ j / ω j , ..., κ N j / ω N j ] ′ and c j = [ C j , ..., C N j ] ′ , with variance matrix Ω j = dia g ( ω j ) . Ifwe elicit a Gaussian prior distribution for the slope coefficients, with p ( β j ) = N (cid:18) µ β j , Σ β j (cid:19) , theconditional posteriors for the slope coefficients are also Gaussian: p (cid:16) β j | β − j , ρ, ω j (cid:17) = N (cid:16) µ β j , Σ β j (cid:17) (3.3) µ β j = Σ β j (cid:20) (cid:16) A − j X (cid:17) ′ (cid:0) κ j − Ω j c j (cid:1) + Σ β j − µ β j (cid:21) Σ β j = (cid:16) A − j X (cid:17) ′ Ω (cid:16) A − j X (cid:17) + Σ β j − . (3.4)The Gaussian conditional posterior of the slope parameters reveals the particular appeal of usinglatent Pólya-Gamma distributed variables. A wide variety of Bayesian model extension, such asvariable selection, or uncertainty over the W can be easily introduced in the above framework.Following Polson et al. (2013), the conditional distribution of ω j is also a Pólya-Gammadistribution: p (cid:0) ω j | β , ..., β J , ρ , ..., ρ J , ω − j (cid:1) = PG (cid:16) , η j (cid:17) , (3.5)where η j = [ η , ..., η N ] ′ . Computationally efficient algorithms for sampling from the Pólya-Gammadistribution are readily available in the R package BayesLogit .The conditional posterior of ρ relates directly to the multinomial logit: p (cid:0) ρ j | ω , ..., ω J , β , ..., β J (cid:1) ∝ p ( ρ j ) N Ö i = J Ö j = (cid:0) exp µ ij (cid:1) y i j Í Jj ′ = exp µ ij ′ (3.6) µ j = A − X β j (3.7)where p ( ρ j ) denotes the prior distribution of ρ j . The conditional posterior in Eq. (3.6) is not froma well-known form and thus cannot be sampled from easily. This is usual in the spatial econometricliterature, and the standard solution is to use a Metropolis-Hastings step, as in LeSage and Pace(2009). 6 arkov-chain Monte Carlo sampling procedure Given the conditional posterior distributions stated above, Markov-chain Monte Carlo algorithmscan be employed by sequentially sampling from the conditional posteriors. We follow the usualidentification assumption of the multinomial logit model in that we set β J = ρ J =
0, and ω J =
0. With suitable starting values for β , ..., β J − and ρ , ..., ρ J − , our sampler involves thefollowing steps:I. For j = , ..., J −
1, update ω j by drawing from p (cid:0) ω j | β , ..., β J , ρ , ω − j (cid:1) using Eq. (3.5),II. For j = , ..., J − β j by drawing from p (cid:16) β j | β − j , ρ , ω j (cid:17) using Eq. (3.3),III. Update ρ j using a Metropolis-Hastings step from p (cid:16) ρ j | ω , ..., ω J , ρ − j , β , ..., β J (cid:17) based onEq. (3.6).The Markov-chain Monte Carlo algorithm cycles through steps I to III B times by excluding thefirst B draws as burn-ins. Inference on the parameters is conducted using the B − B remainingdraws. In a Monte Carlo study we benchmark the SAR multinomial logit model in order to assess thepredictive performance of our proposed modelling framework against two competing specifications:(i) a non-spatial version of the SAR multinomial logit, where all spatial autoregressive coefficients ρ j = j , and (ii) J − J . For the simulation study we use a SAR multinomial logit model as a benchmark data generatingprocess, with three choice classes ( J =
3) and two randomly generated explanatory variables( k = y ij = exp ˜ µ ij Í Jj ′ = ˜ µ ij ′ (4.1)where ˜ µ j = (cid:0) I − ρ j ˜ W (cid:1) − ˜ X ˜ β j ˜ β − J = . . ! + N , − . − .
25 1 ! ! , ˜ β − J = [ ˜ β , ˜ β ] , and ˜ β J =
0. The slope coefficients and the explanatory variables are generatedanew in each Monte Carlo iteration, where ˜ X stems from a standard normal distribution. The Convergence of the MCMC algorithm was checked using the convergence diagnostics proposed by Geweke (1992)and Raftery and Lewis (1992). Convergence diagnostics have been calculated using the R package coda . The SAR logit model is estimated using the method put forward in Krisztin and Piribauer (2020). W isbased on a random spatial pattern generated from a Gaussian distribution for latitude and longitude,and constructed using seven nearest neighbours. Note, that our dependent variable ˜ y ij is a sharevariable, as is often used in land use share models (see, e.g. Chakir and Parent, 2009).To assess the strength of the specifications along multiple scenarios, we vary the strength ofspatial dependence ρ j ∈ { , . , . } . To evaluate the accuracy of the sampler with respect to thechosen sample size, we consider N ∈ { , } . Across all models, our prior set up is as follows:we use a rather uninformative Gaussian prior for β , ..., β J − with zero mean and variance 10 and for ρ , ..., ρ J − we use a the standard beta prior specification as proposed in LeSage and Pace(2009). Table 1:
Root mean squared error measures for the Monte Carlo runs N Model RMSE ρ j = . ρ j = . ρ j = . ρ j direct indirect ρ j direct indirect ρ j
400 SAR Multinomial logit 0.018 0.024 0.192
Multinomial logit
Multinomial logit 0.011
Notes:
Results are based on 1,000 Monte Carlo runs. For each Monte Carlo run, the corresponding sampling algorithms are runusing 1,000 draws, where the initial 700 draws were discarded as burn-in. The columns direct and indirect correspond to summarymarginal effects (for details, see the Appendix). The values given for direct , indirect , and ρ corresponds to the average RMSE (·) over all Monte Carlo iterations. Bold values denote the lowest average RMSE scores. The results of the Monte Carlo study are summarized in Table 1. Each element of the tablecorresponds to the average over 1,000 runs for a particular model specification and Monte Carloscenario. The first and second columns contain information on the sample size N and the modelspecifications. Corresponding to the choice of spatial dependence, the table reports the averageroot mean squared error (RMSE) point estimates for average direct and indirect impacts, as wellas average estimates for ρ j for all j .In the case of no spatial autocorrelation ( ρ j = N = , ρ j = . N = ρ j = . Recent literature focused attention to land sealing resulting from urban sprawl, and associatedspillovers to other land use classes. Results from van Vliet (2019) suggest that in the last decadein Europe 8.4 Mha of land has been converted to urban, out of which 6.3 Mha was converted fromcropland. However, this land sealing led to 13.1 Mha displacement of other land use classes, ascropland was expanded elsewhere, to compensate for the lack of production resources, out of whichthe majority (13 Mha) was expanded in other regions. These spillover effects are well documentedin the literature (e.g. Coisnon et al., 2014; Guastella et al., 2017; Zoppi and Lai, 2014), and serveas a motivation for an empirical application of the spatial multinomial logit model. Both global(Ay et al., 2017) and local (Deng et al., 2008) spillovers are considered of importance.In the spirit of Chakir and Parent (2009), Zoppi and Lai (2014), and Lai and Zoppi (2017) wemodel the areal share of urban sprawl stemming from non-urban land in a given region within aspatial Durbin multinomial logit model, where the log odds take the following form: µ j = ρ j W µ j + α + X β j + W X θ j + ε j . (5.1)The scalar α is an intercept and the term W X is a spatial lag of the matrix of covariates withassociated vector of parameters θ j . This lag explicitly controls for the regions’ characteristics oftheir neighbors. Our sample covers a cross-section of 1,316 European regions across 27 countries. The regionsare classified under the NUTS 2013 classification at the NUTS 3 level. They vary in size andpopulation, however, they divide the territory of the EU for the purpose of harmonized regionalstatistics and analysis. Further, they are assumed to be appropriate spatial observation units foreconomic research and regional policy applications. The regions included in the sample are locatedin Austria (35 regions), Belgium (44 regions), Bulgaria (28 regions), Cyprus (one region), CzechRepublic (14 regions), Denmark (eleven regions), Estonia (five regions), Finland (19 regions),9rance (96 regions), Germany (402 regions), Greece (52 regions), Hungary (20 regions), Italy(110 regions), Latvia (six regions), Lithuania (ten regions), Luxembourg (one region), Malta (tworegions), Netherlands (40 regions), Poland (72 regions), Portugal (25 regions), Republic of Ireland(eight regions), Romania (42 regions), Slovakia (eight regions), Slovenia (twelve regions), Spain(59 regions), Sweden (21 regions), and the United Kingdom (173 regions).The dependent variable for our analysis describes the areal share of urban sprawl emanatingfrom any non-urban type of land within the period from 2000 to 2018. More formally, it is defined asthe land area of a certain type of land use that is being transformed to urban land use between 2000and 2018, divided by the whole area of urban expansion that took place in the respective period.As a result, we obtain a compositional data vector that – by definition – sums up to unity. The typesof land use we consider follow the empirical literature on land use changes and urban expansion(Chakir and Parent, 2009; Chakir and Le Gallo, 2013; Lai and Lombardini, 2016; Zoppi and Lai,2014; Lai and Zoppi, 2017). We distinguish between the five classes cropland, grassland, forest,other, and urban. At this point it is worth noting that we recognize all artificial surfaces as urbanregion, as we want to focus our study especially on soil sealing. The raw data stem from theCORINE Land Cover (CLC) maps provided by Copernicus Land Monitoring Service (CLMS).Their maps are based on satellite data with Minimum Mapping Units (MMU) of 25 hectares forareal phenomena and a minimum width of 100 meter for linear phenomena. The data consists ofan inventory of land cover in 44 classes, which we summarize to the five classes stated above. We use CLC change-layers also provided by CLMS, designed to capture the land cover changesat a higher resolution between two neighbour surveys. Regional aggregates at the NUTS 3 levelare obtained by simple summation of all changes of the corresponding raster elements. Likewise,changes for the whole investigated period are obtained by addition of the 3 sub-periods for whichCLC change-layers are provided. Further data sources are i) the Urban Data Platform Plus providedas a joint initiative of the Joint Research Centre (JRC) and the Directorate General for Regionaland Urban Policy (DG REGIO) of the European Commission, ii) Eurostat (the statistical office ofthe European Union) and iii) the European Observation Network for Territoral Development andCohesion (ESPON).Our set of covariates consists of K ′ =
19 candidate variables that are commonly employedin the literature on land use changes (for an overview, see Shaw et al., 2020). Further, to capturethe complex spatial structure we include not only the spatially lagged dependent vector, but alsothe spatially lagged forms of the explanatory variables (except for the dummy variables). We alsoinclude a vector of ones as intercept. Therefore, the resulting design matrix is of column-dimension K = K ′ + + =
38 where 18 are the spatially lagged covariates and 1 is the intercept. Table 2provides a short technical description for the variables included in our estimation.Since the rent of a certain land use class is assumed to affect the decision of land-owners –yet it is usually not observed – many recent studies consider various proxies to control for thevariation in returns from different land uses (see, e.g. Livanis et al., 2006; Lubowski et al., 2008).Chakir and Parent (2009) conclude that agricultural gross value added divided by the respective Table A1 in the appendix summarizes how each of the land cover classes was mapped. able 2: Variables used in the empirical analysis
Variable Description Source
Croplandto artificial Sum of 2000-2006, 2006-2012, and 2012-2018 CLC land-coverchanges from cropland to artificial land, divided bythe total change of artificial area in the same period. CLCForestto artificial Sum of 2000-2006, 2006-2012, and 2012-2018 CLC land-coverchanges from forests to artificial land, divided bythe total change of artificial area in the same period. CLCGrasslandto artificial Sum of 2000-2006, 2006-2012, and 2012-2018 CLC land-coverchanges from pastures and grassland to artificial land, dividedby the total change of artificial area in the same period. CLCOtherto artificial Sum of 2000-2006, 2006-2012, and 2012-2018 CLC land-coverchanges from area of other use to artificial land, divided bythe total change of artificial area in the same period. CLCCrop rent Share of agricultural gross value added, divided by square kmof area used to grow crops, 2000. JRC, CLCForest rent Share of agricultural gross value added, divided by square kmof forest-area, 2000. JRC, CLCGrass rent Share of agricultural gross value added, divided by square kmof pasture and grassland, 2000. JRC, CLCInitial artificial area Area of artificial land cover, 2000. CLCArtificial growth Growth of artificial areas between 2000 and 2018measured in percent. CLCEmployment primary Share of employment in the primary sector (NACE A), in totalemployment, 2000. JRCEmployment tertiary Share of employment in the tertiary sector (NACE F to Q) intotal employment, 2000. JRCGdp per capita Gross domestic product divided by population, 2000. JRCPopulation density Population per square km, 2000. JRCElevation Average elevation in meters. CopernicusSlope Average slope in degree. CopernicusSoil moisture Content of liquid water in a surface soil layer of 2 to 5 cm depthexpressed as qubic m water per qubic m of soil, 2000. CopernicusN2000 cropland Share of protected area used to grow crops over total area usedto grow crops, 2000. Natura 2000N2000 forest Share of protected area of forests over total area of forests, 2000. Natura 2000N2000 grassland Share of protected area of pastures and grassland over total areaof pastures and grassland, 2000. Natura 2000N2000 other Share of protected area of other use over total areaof other use, 2000. Natura 2000Objective 2 region Dummy varible, 1 denotes region eligible under objective 22000–2006, 0 otherwise. ESPONFarm density Number of farms divided square km, measuredon NUTS2 level, 2000. EurostatFarm size Total farm area devided by number of farms, measuredon NUTS2 level, 2000. Eurostat W we suppose a neighbourhood structure known as sevennearest neighbour specification, where every region is constrained to be a neighbour of its sevenclosest regions. Our results, however, prove robust to variations in the assumed spatial dependencestructure. 12 .2 Empirical results This subsection presents the Markov Chain Monte Carlo (MCMC) results obtained from 10,000posterior draws for our spatial multinomial logit specification, where the first 5,000 were discardedas burn-in. Straightforward interpretation of coefficient estimates in spatial models could lead todeceptive or misleading conclusions (see, e.g. Anselin, 1988; LeSage and Fischer, 2009). Onepossibility is to provide summary metrics in form of direct, indirect (spillover) and total effects.Following (LeSage and Pace, 2009) we present marginal effects in Table 3. Direct effects are thento be interpreted similar to regular slope coefficients. In turn, indirect effects account for theimpacts due to changes in other regions and are therefore to be interpreted as spillover effects. Wefind a significant class-specific spatial parameter ρ j for cropland as well as grassland, highlightingthe necessity of incorporating the spatial dependence structure in the model. This result confirmsthe findings of Guastella et al. (2017) and especially of van Vliet (2019), in that the expansion ofartificial surfaces on productive land leads to further spillover land conversions in surroundingregions.In addition, the table reports the McFadden pseudo R , which serves as a measure of thegoodness of fit in limited dependent variable models. McFadden (1974) highlights that valuesbetween 0.2 and 0.4 already indicate a rather good fit. The rest of the reported results is to beinterpreted as follows: For a certain explanatory variable the four values reported in the direct effectcolumn represent the class specific (cropland, forest, grassland, and other) responses to variationin the respective explanatory. These responses are the changes of the probabilities to convert therespective class in that region to artificial area. Similar, the four class specific columns for theindirect effect are the responses to changes in the explanatory variable in all other regions.The direct effects of the three types of land rent proxies (crop, forest and grass) confirm resultsfrom Chakir and Lungarska (2017) and Chakir and Parent (2009), in that for each land use classhigher rents imply a significantly lower chance of conversion. Additionally, the joint modelling ina multinomial model indicates that significant spillover effects to other classes are present. Mostnotably, an increase in cropland rents in a region, would also increase the conversion of grasslandto artificial areas. We find analogous relationships for forest rent and cropland, as well as grassrent and cropland.A higher initial level of artificial areas indicates that land sealing of other natural vegetationin the own region has a significantly higher probability as compared to land sealing of the otherland covers under scrutiny. Burnett (2012) have similar findings, where urbanization is a processwhich enforces itself. Moreover, as the crop, grass, and forest land surrounding cities is frequentlythe most productive (Shaw et al., 2020), it seems intuitive that urban expansion would take fromthe comparatively less productive other natural vegetation.The change of artificial area appears tohave no direct or indirect impacts on the allocation of its origin.Regarding the sectoral mix of employment, our results indicate that a higher share of ter-tiary employment in the own region implies a significantly higher probability of other naturalvegetation being converted to artificial land. This reflects the findings of Salvati (2016) and Convergence of the sampler was checked using the diagnostics by (Geweke, 1992) able 3: Summary impact measures for artificial area expansion from each land use class
Direct IndirectCropland Forest Grass Other Cropland Forest Grass OtherCrop rent -0.046 -0.003 -0.055 0.042 -0.002 0.012Forest rent -0.008 -0.014 0.021 -0.038 0.008 0.010Grass rent -0.011 -0.031 -0.006 -0.023 0.004 0.008 0.012Initial artificial area -0.041 -0.003 0.008 -0.019Employment tertiary -0.048 0.002 0.013 -0.056 -0.019 -0.005Gdp per capita 0.022 0.024 -0.065 -0.034 0.046 -0.050 -0.039 -0.015 0.037 0.021Elevation 0.006 -0.017 0.020 -0.005 0.014 -0.006 0.006 -0.017Slope -0.042 0.021 -0.003 -0.002 0.012 -0.051 0.034Soil moisture -0.012 0.009 0.010 -0.005 -0.043 -0.004
N2000 cropland -0.015 0.009 0.007 -0.001 -0.062 0.007 0.021 0.036N2000 forest -0.019 -0.007 -0.011 -0.059 -0.029 -0.052
Objective 2 region -0.108 -0.040 ρ j R Notes : Summary metrics are based on 10,000 Markov Chain Monte Carlo iterations, where thefirst 5,000 were discarded as burn-in. Bold written estimates indicate statistical significance undera 90% credible interval.Salvati and Carlucci (2016), where higher tertiary employment is found to mainly reflect thepresence of urban fabric. In this context, the positive spillover effects of primary and tertiaryemployment to neighbouring regions’ grassland can be contextualized as the effect of industrialbelts on pastures. This reflects the findings of the theoretical model of Turner (2005), where chieflyindustry clusters agglomerate with large-scale livestock farms.Our results with regards to gross domestic product per capita suggest that it is not a significantdriver of land sealing in a European context. This is as opposed to findings of e.g. Deng et al.(2008) in developing countries, where gdp per capita is found to be one of the main driversof urbanization. Moreover, we find that a higher gdp per capita in fact significantly lowers theprobability of sealing grass land in the own region. When observed jointly with the direct effectsof population density, this supports findings by McGrath (2005) and Guiling et al. (2009), whofind that population is a more significant driver of urbanization, as opposed to personal income.14ote, however that only for the land take from grass land is the effect positive and significant. Forland takes from forest and other natural land, a higher population density in fact results in a lowerchance of land conversion. This result can be interpreted on the one hand with the fact that regionswith a higher endowment of population density are more urban in nature and contain a much lowerpercentage of cropland or other natural vegetation. On the other hand, specific literature, such asDelbecq et al. (2014) and Wu (2006), provide evidence that private home-owners exhibit strongpreferences for surrounding grassland amenities.Turning our attention to the estimated impacts of the biophysical drivers elevation, slope,and soil moisture, we can largely confirm the overall conclusions of Shaw et al. (2020) andChang-Martínez et al. (2015) in that the biophysical processes play a secondary role to socio-economic ones in explaining land sealing processes. For the own-region, only slope plays has asmall, albeit significant impact on the probability of sealing other natural vegetation. Additionally,a higher percentage of soil moisture indicates a significantly higher chance of converting grasslandto urban land in neighbouring regions.Our results seem to indicate that the
Natura 2000 protection program has intended effects,as higher shares of protected forest and grassland would – according to our results – lead tosignificantly lower chance of the respective land cover being converted into urban. Note, that if aregion has a higher share of other natural vegetation under Natura 2000 protection, this would lowerthe chances of neighbouring regions converting this land cover to urban. This largely confirms thefindings of Lai and Lombardini (2016), Zoppi and Lai (2014), and Lai and Zoppi (2017). Our jointmultinomial logit framework, however, allows us to uncover additional interdependencies amongthe natural protection of land covers. Our estimated results suggest that a higher share of protectedforest in a NUTS3 level region would results in a significantly higher probability of convertingcropland to urban, not only in the own region but also amongst neighbours. Additionally, thiswould also lower the odds sealing grassland under artificial surfaces.The estimated results with regard to our subsidy proxies seem to show that regional fundingplays a comparatively larger role as farm-specific subsidies. The own-regional effect of regionallevel Objective 2 subsidies is significant and negative for cropland, and positive for other nat-ural vegetation. This finding supports the hypotheses that subsidies increase land conversion(Shaw et al., 2020). Particularly noteworthy is the result that the land take comes more signifi-cantly from natural vegetation (which is highest in biodiversity) as opposed to more productivecropland. Additionally, neighbours of regions under Objective 2 funding also have a significantlydecreased chance of converting cropland to artificial areas. This might suggest an increase incropland productivity, through better infrastructure. With regard to our farm structure variables –proxying the role of the Common Agricultural Policy – our results indicate a significant effect onlywith regard to farm density: a higher density of farms would lower the chance of converting landto grassland in the own region. 15
Concluding remarks
In this paper we put forth a Bayesian estimation approach for a multinomial logit specification forthe modelling of land use conversion, which has a spatial autoregressive structure in the log odds,with differing strength of spatial autocorrelation for each choice alternative. The virtue of ourspecification is that it combines a spatial autoregressive framework (allowing for cross regionalspillovers), and a joint multinomial framework (allowing for cross land use class dependencies).The proposed approach is based on recent spatial econometric advances dealing with Bayesianestimation of the logit model Krisztin and Piribauer (2020). The core step of the estimationprocedure relies on introducing a latent Pólya-Gamma variable (see Polson et al., 2013). Throughthe latent variable, the conditional posterior distribution of the slope parameters in the spatialautoregressive logit specification is rendered in a Gaussian form, which allows us to tackle theMCMC estimation in a particularly efficient way. We demonstrate in a simulation study theadvantages and behaviour of our proposed model specification, benchmarking it against simpleralternatives.The virtues of the spatial multinomial logit model are illustrated using an empirical specifica-tion. Specifically, we examine the land sealing activities in European NUTS-3 level regions. Weconsider the areal share of urban sprawl emanating from cropland, grassland, forest, and other nat-ural vegetation from 2000 to 2018. The observation on land use data stem from the CORINE LandCover (CLC) maps. Our results suggest, that spatial dependence indeed play a small, but signifi-cant role, particularly for the land use classes cropland and grassland. For all land covers proxiedland rents are of central importance. Additionally, our findings corroborate evidence from recentliterature, that socio-economic drivers play a much more central role, as opposed to biophysicalones (for an overview, see Shaw et al., 2020). The key role of population density Guastella et al.(2017); Deng et al. (2008); Lai and Lombardini (2016) in urban land take is confirmed by ourresults. Moreover, we confirm on a larger level that environmental protection not only has effectsin the own- but also in neighbouring regions (Lai and Lombardini, 2016; Zoppi and Lai, 2014;Lai and Zoppi, 2017). Through the virtue of our multinomial analysis we also find evidence forforest protection having spillover effects to neighbouring regions and other land covers.
References
Alexiadis S, Hasanagas N and Christos L (2013) Examining Agriculture from a Regional Perspec-tive: Implications for the Common Agricultural Policy.
Land Use Policy
Spatial Econometrics: Methods and Models . Kluwer Academic, DordrechtArtmann M (2014) Assessment of soil sealing management responses, strategies, and targetstoward ecologically sustainable Urban land use management.
Ambio
Environmental Modeling & Assessment
Land Economics
Land Economics
Ecological Economics
Spatial Economic Analysis
Papers in Regional Science
International Journal ofGeo-Information
Regional Science and Urban Economics
Land Economics
Journal of Urban Economics
Austrian Journal of Statistic
Statistics and Computing
BayesianStatistics . Oxford University Press, Oxford, 4th edition, 167–193Guastella G, Pareglio S and Sckokai P (2017) A spatial econometric analysis of land use efficiencyin large and small municipalities.
Land Use Policy
63, 288–297Guiling P, Brorsen BW and Doye D (2009) Effect of urban proximity on agricultural land values.
Land Economics
Bayesian Analysis
American Journal of Agricultural Economics
Global Assessment Report on Biodiversity and Ecosystem Services . IPBES, BonnKlier T and McMillen DP (2008) Clustering of auto supplier plants in the United States: Generalizedmethod of moments spatial logit for large samples.
Journal of Business and Economic Statistics
Empirical Economics
Lai S and Lombardini G (2016) Regional drivers of land take: A comparative analysis in two17talian regions.
Land Use Policy
56, 262–273Lai S and Zoppi C (2017) The influence of natura 2000 sites on land-taking processes at theregional level: An empirical analysis concerning Sardinia (Italy).
Sustainability (Switzerland)
Preprint
LeSage JP and Fischer MM (2009) Spatial growth regressions: Model specification, estimationand interpretation.
Spatial Economic Analysis
Journal of the Royal Statistical Society. Series A:Statistics in Society
Introduction to Spatial Econometrics . CRC Press, Boca RatonLondon New YorkLi M, Wu JJ and Deng X (2013) Identifying drivers of land use change in China: A spatialmultinomial logit model analysis.
Land Economics
American Journal of Agricultural Economics
LandEconomics
Frontiers of Econometrics . Academic Press, New York, 105–142McGrath DT (2005) More evidence on the spatial scale of cities.
Journal of Urban Economics
AmericanJournal of Agricultural Economics
Regional Science and Urban Economics
Journal of the American Statistical Association
Bayesian Statistics . Oxford University Press,Oxford, 4. edition, 763–773Salvati L (2016) Soil sealing, population structure and the socioeconomic context: a local-scaleassessment.
GeoJournal
Urban Studies
Landscape and Urban Planning
Regional Environmental Change
Journal ofUrban Economics
Nature Sustain-ability
Journal ofEnvironmental Economics and Management
Journal of Environmental Economics and Management
Land Use Policy ppendixMarginal effects
Similar to the marginal effects of the spatial Durbin logit model (see Krisztin and Piribauer, 2020),in the presented multinomial logit model in Eq. (5.1) the interpretation of marginal effects of the k -th explanatory variable (with k = , ..., K ) differs from those in linear models. This is due to thefact that the multinomial logit model is non-linear in nature, but also the the presence of spatialautocorrelation gives rise to an N × N matrix of partial derivatives, which makes interpretation ofmarginal effects richer, but also more complicated (see also LeSage and Pace 2009).As is standard in the logit literature, and analogous with the proposed marginal effects of thespatial logit model (Krisztin and Piribauer, 2020), we provide marginal effects relative to the meanof the k -th explanatory variable, which we denote as x k = Í Ni = x ik / N . Thus the interpretation ofthe marginal effects is the change in probability of observing y = j associated with a change in theaverage sample observation of the k -th explanatory variable. To write the partial derivatives of themodel in Eq. (5.1), with respect to the k -th coefficient let us define: µ k j = A − j I N x k β k j + A − j W x W k θ k j , ζ k j = A j − I N β k j + A − j W θ k j , and p k j = exp µ k j Í Jj ′ exp (cid:16) µ k j ′ (cid:17) .β k j and θ k j denote the k -th element of β j and θ j , respectively. x W k denotes the average value ofthe k -th spatially lagged explanatory variable. The partial derivatives can then be expressed as: ∂ p ( y = j | x k ) ∂ x ′ k = p k j ⊙ " ζ k j − J Õ j ′ p k j ′ ⊙ ζ k j ′ , (A.1) = Λ k j , where ⊙ is the Hadamard product. Note that marginal effects of the k -th coefficient on class j ,denoted as Λ k j , are an N × N matrix due to the presence of the N × N spatial multiplier A − j .Interpreting N × N marginal effects proves cumbersome, therefore we define summary impacteffects (LeSage and Pace, 2009). These can be readily calculated from Λ k j : direct k j = N ι ′ N dia g ( Λ k j ) (A.2) total k j = N ι ′ N Λ k j ι N (A.3) indirect k j = total k j − direct k j , (A.4)where ι N denotes an N × able A1: Mapping of detailed CLC classes to land use aggregates
CLC code Description Aggregate
111 Continuous urban fabric Artificial112 Discontinuous urban fabric Artificial121 Industrial or commercial units Artificial122 Road and rail networks and associated land Artificial123 Port areas Artificial124 Airports Artificial131 Mineral extraction sites Artificial132 Dump sites Artificial133 Construction sites Artificial141 Green urban areas Artificial142 Sport and leisure facilities Artificial211 Non-irrigated arable land Cropland212 Permanently irrigated land Cropland213 Rice fields Cropland221 Vineyards Cropland222 Fruit trees and berry plantations Cropland223 Olive groves Cropland231 Pastures Grassland241 Annual crops associated with permanent crops Cropland242 Complex cultivation patterns Cropland243 Land principally occupied by agriculture & natural vegetation Forest244 Agro-forestry areas Forest311 Broad-leaved forest Forest312 Coniferous forest Forest313 Mixed forest Forest321 Natural grasslands Other322 Moors and heathland Other323 Sclerophyllous vegetation Other324 Transitional woodland-shrub Other331 Beaches dunes sands Other332 Bare rocks Other333 Sparsely vegetated areas Other334 Burnt areas Other335 Glaciers and perpetual snow excluded415 Inland marshes excluded412 Peat bogs excluded421 Salt marshes excluded422 Salines excluded423 Intertidal flats excluded511 Water courses excluded512 Water bodies excluded521 Coastal lagoons excluded522 Estuaries excluded523 Sea and ocean excluded999 NODATA excluded111 Continuous urban fabric Artificial112 Discontinuous urban fabric Artificial121 Industrial or commercial units Artificial122 Road and rail networks and associated land Artificial123 Port areas Artificial124 Airports Artificial131 Mineral extraction sites Artificial132 Dump sites Artificial133 Construction sites Artificial141 Green urban areas Artificial142 Sport and leisure facilities Artificial211 Non-irrigated arable land Cropland212 Permanently irrigated land Cropland213 Rice fields Cropland221 Vineyards Cropland222 Fruit trees and berry plantations Cropland223 Olive groves Cropland231 Pastures Grassland241 Annual crops associated with permanent crops Cropland242 Complex cultivation patterns Cropland243 Land principally occupied by agriculture & natural vegetation Forest244 Agro-forestry areas Forest311 Broad-leaved forest Forest312 Coniferous forest Forest313 Mixed forest Forest321 Natural grasslands Other322 Moors and heathland Other323 Sclerophyllous vegetation Other324 Transitional woodland-shrub Other331 Beaches dunes sands Other332 Bare rocks Other333 Sparsely vegetated areas Other334 Burnt areas Other335 Glaciers and perpetual snow excluded415 Inland marshes excluded412 Peat bogs excluded421 Salt marshes excluded422 Salines excluded423 Intertidal flats excluded511 Water courses excluded512 Water bodies excluded521 Coastal lagoons excluded522 Estuaries excluded523 Sea and ocean excluded999 NODATA excluded