Large scale spatio-temporal patterns of road development in the Amazon rainforest

Summary

There is burgeoning interest in predicting road development because of the wide ranging important socioeconomic and environmental issues that roads present, including the close links between road development, deforestation and biodiversity loss. This is especially the case in developing nations, which are high in natural resources, where road development is rapid and often not centrally managed. Characterization of large scale spatio-temporal patterns in road network development has been greatly overlooked to date. This paper examines the spatio-temporal dynamics of road density across the Brazilian Amazon and assesses the relative contributions of local versus neighbourhood effects for temporal changes in road density at regional scales. To achieve this, a combination of statistical analyses and model-data fusion techniques inspired by studies of spatio-temporal dynamics of populations in ecology and epidemiology were used. The emergent development may be approximated by local growth that is logistic through time and directional dispersal. The current rates and dominant direction of development may be inferred, by assuming that roads develop at a rate of 55 km per year. Large areas of the Amazon will be subject to extensive anthropogenic change should the observed patterns of road development continue.

INTRODUCTION

Roads are an important part of everyday life for most people. Road development influences a wide range of phenomena, from human society, business and economies, to the natural environment (Forman et al. 2003). In regional development, roads are often perceived as the initial stage of development, especially in the tropics where they open access to remote areas for colonization, agricultural extraction and resource extraction (Laurance et al. 2001; Arima et al. 2005; Perz et al. 2007; Caldas et al. 2010). In contrast to the predominantly positive influences roads have on human enterprises, they have many and varied effects on the environment. Examples include: fragmenting habitats and altering their structure, altering animal behaviour, movement patterns and habitat use, and introducing pollutants (for more details see Forman 1998; Forman & Alexander 1998; Spellerberg 2002; Forman et al. 2003; Coffin 2007; Laurance et al. 2009). The majority of these impacts are negative, with one estimate suggesting that impacts on fauna are five times more likely to be negative than positive (Fahrig & Rytwinski 2009) and Benitez-Lopez et al. (2010) showed that the impact is greater closer to roads. In the Amazon, much work has been done to study the link between roads and deforestation, with findings highlighting the close spatial link between road development and deforestation rates (Geist & Lambin 2002; Walker et al. 2004; Perz et al. 2007; Fearnside 2008). Banning road development in important and delicate tropical areas has been suggested as a way to prevent these negative impacts (Laurance et al. 2009). However, this is highly unfeasible given the socioeconomic benefits and development potential that roads bring (Maki et al. 2001). As a result, Laurance and Balmford (2013) suggested global road zoning exercises to determine the future location and design of road networks, aiding policy makers to maximize benefits to people while minimizing costs to the environment. Ideally, such exercises should not only incorporate existing information on road networks, but also information on how road development changes through space and time in response to different natural and anthropogenic factors. Our study attempts to understand and predict patterns of road development at large spatial scales, characterizing large-scale spatio-temporal patterns of road development in the Brazilian Amazon using predictive models of road development.

One complication in understanding the patterns of road development is that roads are constructed for many and varied purposes. As a result, several distinctions between different road types have aided in understanding their characteristics. Examples include: exogenous versus endogenous, where endogenous roads are built as a result of economic development and exogenous roads are built to stimulate economic development; paved versus unpaved, which distinguishes roads based on their physical characteristics; primary versus secondary, with primary roads referring to large highways and secondary roads referring to roads that build off from primary roads; and official versus unofficial roads, which distinguishes roads based on whether planning application for the construction of the roads was provided (Perz et al. 2007; Pfaff et al. 2007). Official roads are built either by the government or with government permission, and tend to be large-scale projects such as the Transamazonian highway, which have planning put into place years in advance of any construction. Unofficial roads are built with no planning permission obtained from the state, often by non-state actors, such as miners or loggers, and tend to be built on a smaller scale than official roads (Brandão & Souza 2006; Perz et al. 2007). In the Brazilian Amazon, most primary roads were built in the 1970s, and they tend to be official. An extensive network of secondary roads (mostly unofficial) has developed over the years; because the majority of roads built are unofficial there is a distinct lack of spatial information on their location and extent (Brandão & Souza 2006). This presents a problem for policy makers and conservationists who need spatial information on current and future roads in order to assess potential impacts and make informed decisions.We focus on characterizing the spatio-temporal dynamics of secondaryroads at an Amazon-wide scale.

Predicting and understanding road development is made particularly challenging because of the many factors that influence road development. At large scales, economics, policy, technology, demographic and cultural factors, influence the rate, location and extent of road development (Geist & Lambin 2002; Montagnini & Jordan 2005). For example, government policies greatly influence investment in roads; government may provide subsidies to road builders or directly invest in the network by building federal roads to stimulate economic growth. This process was exemplified by the drive to colonize new areas in the Brazilian Amazon in the 1970s through building long-distance official roads (Carvalho et al. 2002). At that point in time the motivation for building roads was to connect existing cities and secure frontier regions, little thought was given to constraining factors such as topography. At smaller scales, road alignment (the location of the road in relation to the surroundings) is dependent on a range of factors that can either constrain or facilitate the laying of the road, including topography, existing land use, hydrological features and ground conditions (Koorey 2009).

Despite the challenges in predicting road development, it remains a major topic of research because of the myriad socioeconomic and environmental consequences (Arima et al. 2005). For example, road development predictions are used to aid environmental impact assessments, producing likely scenarios of environmental change and associated impact estimations, such as predicting future biodiversity outcomes (Soares-Filho et al. 2006; Alkemade et al. 2009). However, as described, multiple factors influence road development, and potentially do so in complicated ways involving nonlinear feedbacks and interactions, time lags and stochasticity (Ralston & Barber 1982). This resulted in early road models that were focused on describing and explaining the spatio-temporal development of roads, highlighting the influence of geopolitics on road location (Andersen & Reis 1997; Carvalho et al. 2002; Perz et al. 2007), but lacking analytical or predictive elements. Models that went beyond mere description were complicated, and often had poor predictive ability (see for example Taaffe et al. 1963; Rene 1965).

More recent models that predict future road networks are often at the spatial and temporal resolution of individual roads(Soares-Filho et al. 2006; Jiang 2007; Arima et al. 2008), although there seems to be a move towards larger scale (or coarser resolution) networks (see for example Walker et al. 2013). Individual scale road models allow the location of specific roads to be determined; however the validity of these predictions is not well established. Further, given the current predictive modelling as applied to individual road development, the characterization of larger scale spatio-temporal patterns has been overlooked. Models characterizing broad scale dynamics can be easier to generalize and incorporate into other analyses, for instance being incorporated into larger models, integrated assessment models, or being applied to other geographic regions. Approximately two-thirds of papers predicting land-use change in the Amazon region use roads as a predictor of future land use. However the vast majority of land-use change models do not currently attempt to predict the future development of the road network, and treat the dynamic development process as static. Models that can easily characterize road development could thus play a vital role in future land-use modelling.

Given the complex dynamics of road development and its interactions with economics, policy, technology, demographic and cultural factors, it could be asked whether these complicated interactions and driving forces have emergent properties that can be used to predict road development and described using simple models? Recently, Ahmed et al. (2013) showed that the dynamics of road density through time in regions of the Amazon could be characterized as a logistic growth curve, where road density initially grows through time at an approximately exponential rate, but slows as road density approaches a maximum for that region, suggesting that the road development within a region does have general aspects to its behaviour. This raises further questions about other general characteristics of large-scale road dynamics. If roads tend to show logistic dynamics within regions, how then are the regional dynamics related across the Amazon region? A priori, we would expect the spatio-temporal dynamics to reflect the general direction of development in the Amazon region from the densely populated regions in the east towards the Amazon rainforest in the west. However, to date, no-one has characterized these dynamics formally.

Here we characterize the dynamics of road density across the Brazilian Amazon, quantifying the spatio-temporal dynamics in terms of the dominant direction and rate of road development. Further, we assess the relative contributions of intrinsic versus neighbourhood effects for temporal changes in road density at regional scales. Intrinsic effects take into account the development of the roads within any given area. Neighbourhood effects take into account the state of the road network in surrounding areas and determine how this influences the network in any given area.

To achieve our objectives, we use a combination of statistical analyses and model-data fusion techniques inspired by studies of the spatio-temporal dynamics of populations in ecology and epidemiology (Hilborn & Mangel 1997; Ferrari et al. 2008; Haynes et al. 2009). Model-data fusion techniques are computational methods that enable quantitative or qualitative details of mathematical models to be inferred; typically the parameters but it may also be functional forms and even entire model components, by conducting formal repeated comparisons between model predictions and data. Such methods are increasingly popular in ecology for two main reasons: they generate the best predictive model under an assumed model structure, given available data and prior uncertainty, and they enable the objective assessment of parameter probability distributions and the relative likelihood of alternative model formulations given empirical evidence. Our study is inspired by the application of such methods in population ecology and epidemiology because these are also commonly concerned with predicting large spatio-temporal patterns of properties that grow and spread over time, although in ecology and epidemiology the loss of individuals is also important. In addition, given that Ahmed et al. (2013) have already shown roads follow logistic growth dynamics similar to those found among many species’ populations, it is natural to extend that analogy to consider the spatio-temporal dynamics of road development as being similar to spatio-temporal population growth (for example like the invasion of an alien species). It is important to note that while models are often developed to enable deeper understanding of what underpins certain phenomena of interest, our goal was primarily to characterize the large-scale spatio-temporal patterns of road development in the Brazilian Amazon, and to investigate how well simple phenomenological models describe it, although we do extendthis to investigate the importance of explicitly distinguishing between local and neighbourhood effects for one of our model formulations. Future work will hopefully build from these findings to obtain new insights into the dynamic mechanisms.

METHODS

Data and general pattern analyses

We used existing data on road density (km roads per km² land area) for 443 municipalities of the Brazilian Amazon in 2004 and 2007 (Fig. 1). These were obtained from IMAZON (Instituto do Homem e Meio Ambiente da Amazônia, http://www.imazon.org.br/), and derived using manual digitization methods on Landsat TM imagery in ArcGIS 9.3 (see Brandão & Souza 2006 for details) for the purpose of studying the link between the building of new roads and deforestation patterns. All roads, official and unofficial, were included in the analyses; only two complete road maps detailing all roads (2004 and 2007) were available. Given that much of the official (primary) road network (namely highways) was established in the 1970s and has changed little since then, the majority of development over our study time frame has taken place in the secondary road network. For the statistical analyses of spatio-temporal patterns, we analysed the data specifying road density in each municipality in 2004 and 2007. For model training and evaluation, we analysed the data specifying road density in each cell of a 100-km² resolution grid for the same years.

Figure 1 (a) The spatial distribution of the change in road density between 2004 and 2007 in the Amazon legal, divided by municipio. Change in road density is concentrated along the arc-of-deforestation. (b) Histogram of log-change in road density between 2004 and 2007. (c) Directional anisotropy radial plot, displaying the extent to which statistically significant spatial correlation in change in road density extends in different directions. ‘Long’ bands indicate correlation extends to greater distances, indicating that road development is moving at a perpendicular angle (differing grey tones are a visual aid to enable differentiation between bands). (d) Example of road development over a three-year period between 2004 (light grey) and 2007 (dark grey).

Extending the finding of Ahmed et al. (2013), that road development is affected by the existing road density of a given municipality (i.e. intrinsic effects), we carried out a general linear model (GLM) analysis (Crawley 2008) to determine whether road density in neighbouring municipalities (i.e. neighbourhood effects) had an influence on road density growth in a given municipality. We calculated Moran's I for the whole study area to investigate road density anisotropy (spatial correlations) among municipalities (Legendre & Legendre 1998). Further, we used directional anisotropy to determine if the development of road networks was moving in any particular direction. Directional anisotropy analysis was implemented at both Amazon-wide and quadrant scales to see if the directional patterns observed at an Amazon-wide scale held at smaller scales.

Models

We consider a set of simple models to characterize the spatio-temporal patterns of road development across the Brazilian Amazon. These are broadly divided into two different types: those that characterize the patterns as spatio-temporal travelling waves (travelling wave model) and those that predict the changes in road density in individual grid cells by explicitly separating contributions from processes intrinsic to the grid cell from processes occurring in neighbouring grid cells. We fit these latter models with and without assuming the existence of neighbourhood effects to see if accounting for neighbourhood effects influences model performance. Although the travelling wave model does not explicitly account for neighbourhood effects, they are implicitly accounted for in its formulation. All models were calibrated and run on the municipio scale and a 100-km grid; for brevity, only results from the 100-km grid are presented, as there was little difference between the two scales For simplicity, the term ‘municipio’ is used in the model descriptions for any individual municipio or grid cell.

Travelling wave model

The first model was a simple phenomenological description of a travelling wave of road density across Brazilian Amazon, with zero or low road density in front of the wave, the highest road density behind the wave, and the moving ‘wave front’ capturing the change in road density through time and space. The model has the form

(1) \begin{equation} \rho _m \left( t \right) = \frac{{K_w }}{{1 + \exp \{ - R[\gamma \left( {\alpha ,\omega } \right) + \tau c]\} }},\end{equation}

where ρ _m is the road density (km km⁻², namely the length of road km, per area km²) in location (municipio) m at time t (years), K_w is the maximum road density behind the wave front (km km⁻²), R scales the rate of change of road density with space (km km⁻²yr⁻¹), γ(α, ω) is the distance from the midpoint of the location to the midpoint of the travelling wave with coordinates (α) when projected along the angle of movement of the travelling wave relative to north (ω), τ is the number of years that have elapsed since 2004 and c is the speed of travel of the travelling wave (km yr⁻¹). This model was fitted to data on the entire Brazilian Amazon together, and also for individual quadrants of the Brazilian Amazon (see Results), to explore evidence for the wave of road development travelling with different properties in different regions.

No neighbourhood effects models

We analysed two models in which we assumed no neighbourhood effects, to contrast with the neighbourhood effects models. The first model changes in road density through time as a logistic growth process (following Ahmed et al. 2013)

(2) \begin{equation} \frac{{d\rho _m }}{{\it dt}} = \rho _m r\left(1 - \frac{{\rho _m }}{K}\right) = {\it Logistic}\end{equation}

where K is the maximum road density per location and r is the maximum rate of growth of road density through time (km km⁻² yr⁻¹). The second model predicts the change in road density through time as an exponential growth process

(3) \begin{equation} \frac{{d\rho _m }}{{\it dt}} = \rho _m r = {\it Exponential}\end{equation}

where the parameters and assumptions are as defined above. This was investigated simply to contrast a model that assumed density dependence in the growth of road density (Eq. 2), with one that did not.

Neighbourhood effects models

We considered a range of different models in which the growth of road density in each location was a combination of local growth and neighbourhood effects. All models have the general form

(4) \begin{equation} \frac{d\rho _m }{dt} = {\rm local\, growth}\,( {\rm logistic} ) \, + \,{\rm neighbourhood}\,{\rm effects}\quad \end{equation}

All of our preliminary analyses supported the use of the logistic model (Eq. 2) rather than the exponential model (Eq. 3) as the local growth process. Our different neighbourhood effects models thus only differ in their ‘neighbourhood effects’ terms (in Eq. 4). We reasoned that neighbourhood influences will probably vary as a function of the relative difference in road density between neighbouring locations (all those that share borders with it). We hypothesized that the larger the difference in road density between neighbouring locations, the larger the pressure will be on the neighbour with the lower road density to increase in road density, and that this effect may not be reciprocal; so a location with high road density neighbouring one with low road density may not feel any neighbourhood effects, but will be dominated by local growth processes. We considered four different neighbourhood effects formulations (here referred to as NEm1–NEm4)

(5) \renewcommand\theequation{(\arabic{equation}, NEm1)} \begin{eqnarray} &&{\it Neighbourhoodeffects_m}\nonumber\\ &&\quad = D_1 \mathop \sum \limits_{j = 1..n} {\rm exp}(\max ( {\rho _j - \rho _m - \tau ,0} ))\qquad\qquad \end{eqnarray}

(6) \renewcommand\theequation{(\arabic{equation}, NEm2)} \begin{eqnarray} &&{\it Neighbourhoodeffects_m} \nonumber\\ &&\quad= D_2 \mathop \sum \limits_{j = 1..n} {\rm exp}\left(\frac{{\max (\rho _j - \rho _m - \tau ,0)}}{{E_{mj} }}\right)\qquad\quad \qquad\end{eqnarray}

(7) \renewcommand\theequation{(\arabic{equation}, NEm3)} \begin{eqnarray} &&{\it Neighbourhoodeffects_m}\nonumber\\ &&\quad = D_3 \mathop \sum \limits_{j = 1..n} \rho _j {\rm exp}(\max (\rho _j - \rho _m - \tau ,0))\qquad\quad\qquad\end{eqnarray}

(8) \renewcommand\theequation{(\arabic{equation}, NEm4)} \begin{eqnarray} &&{\it Neighbourhoodeffects_m}\nonumber\\ &&\quad = D_4 \mathop \sum \limits_{j = 1..n} \rho _j {\rm exp}\left(\frac{{\max (\rho _j - \rho _m - \tau ,0)}}{{E_{mj} }}\right)\qquad\quad\qquad\end{eqnarray}

where D scales the magnitude of the neighbourhood effect (units differ depending on formulation), j is the identity of one of the n neighbours to location m, E_mj is the Euclidean distance (km) between the centroids of locations m and j, and τ is the threshold difference in road density between a neighbour and the focal location, below which the focal location experiences no neighbourhood effects. Equation (5) predicts the size of the neighbourhood effect as an exponential function of the sum of positive differences in road density between a location and its neighbours. Equation (6) assumes that there is an additional effect based on the relative distance between neighbours. Equation (7) assumes that the size of the neighbourhood effect is an exponential function of the difference in road density between a location and that of its neighbour, and of the absolute road density of that neighbour. Equation (8) also assumes this, and an additional effect from the Euclidean distance between neighbouring sites. We settled on the exponential transformation in Equations (5–8) after comparing models assuming linear, exponential and saturating functional transformations; reflecting different hypotheses about the difference neighbourhood density and the size of the neighbourhood effect. These analyses (Appendix 1, see supplementary material at Journals.cambridge.org/ENC) suggested that the exponential function performed best, but only marginally, and so we used that function throughout the rest of our study.

Maximum likelihood parameter estimation

We inferred the probability distributions of the parameters to the above models given the data using maximum likelihood parameter estimation. We expected that most of the variance in the observed change in road density derived from inherent stochasticity in the road development process rather than any error in taking the measurements of road density, which we were confident were accurate (Brandão & Souza 2006; Ahmed et al. 2013). We therefore expected larger variation amongst sites that showed relatively large changes in road density over time due to the multiplicative nature of road density growth. Preliminary analysis of the data also indicated lognormal variation in the change in road density across all sites. Therefore for all models, except the travelling wave model (for reasons detailed below) we estimated the likelihood of the data given the predictions of a parameterized model using

(9) \begin{eqnarray} &&\log \left( {L\left( {{\it Model}\left( \theta \right){\rm |}X} \right)} \right)\nonumber\\ &&\quad = \mathop \sum \limits_{j = 1..J} \log ({\rm P}(\log ( {\rho _j ( {2007} ) - \rho _j ( {2004} )} )\nonumber\\ &&\quad \approx \aleph (\log ( {\rho _{j,pred} } ),\sigma ^2 ))\end{eqnarray}

where L is the likelihood of the model with vector of parameters θ given a vector of data, X of length J, P is the probability that the log of the observed change in road density between 2004 and 2007 is drawn from a normal distribution, ℵ, with mean centred on the log of the predicted change in road density over that time period, Δρ _{j, pred} , and variance σ² which is estimated alongside the other model parameters. We excluded sites that exhibited no change in road density over the time window because the model formulation uses log road density change and we were unable to log zero values. This could potentially bias our parameter estimates by excluding sites exhibiting no change (especially zeros) that we would predict would experience neighbourhood effects. While we cannot rule out this possibility entirely, the fact that only a low fraction of our sites exhibited no change over the time period (43 out of 474 municipios exhibited no change, and only a fraction of these would potentially experience neighbourhood effects) and that we obtained very similar estimates for the rate of spread of road density for the neighbourhood effects and the travelling wave model (where we did not exclude sites exhibiting no changes) indicated to us that any introduced biases were unlikely to have been large.

For the travelling wave model we used the 2004 and 2007 data to infer the most likely position and parameters of a unidirectional travelling wave

(10) \begin{eqnarray} &&\log \left( {L\left( {{\it Model}( \theta ){\rm |}X} \right)} \right)\nonumber\\ &&\quad = \mathop \sum \limits_{\tau = 0,t = 3} \mathop \sum \limits_{j = 1..J} \log ({\rm P}(\log ( {\rho _j ( {2004 + \varepsilon } )} )\nonumber\\ &&\quad \approx \aleph (\log ( {\rho _{j,pred} } ),\sigma ^2 ))\end{eqnarray}

where ε = 0 corresponds to the 2004 road density data and ε = 3 corresponds to the 2007 road density data. The only thing that changes in Equation (10) between those two time points is the ε parameter, which enabled us to infer the wave speed parameter c. Inference of the travelling wave model was performed using the raw data rather than the change in road density, because it enabled us to infer the speed and direction of the travelling wave. This also meant that we did not have to exclude sites with no change in road density over the time window.

For all models we used Markov Chain Monte Carlo sampling with the Metropolis-Hastings algorithm (Gilks et al. 1996) to perform the parameter estimation, which we implemented using the Filzbach libraries (http://research.microsoft.com/en-us/um/cambridge/groups/science/tools/filzbach/filzbach.htm). All Markov Chains were 100000 iterations in length after a 10000 iteration burn in period, and sampled every 100 iterations (the default setting), which was deemed sufficient to remove autocorrelation after preliminary visual and statistical analyses. We used 10-subset cross validation to assess the sensitivity of our inferred parameter probability distribution to different subsets of the data and to assess the performance of the fitted models against data not used in training (to prevent over-fitting).

We assessed the goodness of fit of the model to the data using five different metrics: the deviance information criterion (DIC; Gelman et al. 2004), the correlation between the model predictions and the evaluation data (CC), the coefficient of determination between the model predictions and the evaluation data (CD), and the mean log likelihood of the training (TL) and evaluation data (EL).

Future projections

The model that explained the highest percentage of the variation in the evaluation data was used to make future projections of road network spread. While the results need to be interpreted with caution due to the limited amount of calibration data, especially the use of a single three-year time window over which changes in road density were observed, they highlight how the models can be used to make projections and provide a first estimate using this approach. Further, we modify the model to incorporate the effect of barriers to development (such as protected areas) and investigate how the projections are affected.

RESULTS

General patterns

The spatial distribution of the change in road density was concentrated along the arc-of-deforestation (Fig. 1 a). The change between 2004 and 2007 was approximately log-normally distributed based on visual inspection of the data (Fig. 1 b), although a Shapiro-Wilks test indicated that the distribution was significantly different from log-normal (p < 0.001). A GLM indicated that the growth in road density is significantly positively related to the initial density of a municipio (p < 0.001, slope = 0.11) and also indicated a significant interaction effect between the initial density of a given municipio and its average neighbourhood road density on the growth of road density (p < 0.001, df = 470, slope = –0.36). However, average neighbour density alone was not significantly related to the change in road density (p = 0.18).

Moran's I indicated that there was significant autocorrelation in road density change between municipalities that were up to 434 km apart. The analysis of spatial anisotropy in the change in road density between 2004 and 2007 implies a number of different directions of the development of the road network across the whole Amazon, with correlation with distance dropping away most sharply in both north-westerly and north-easterly directions (between 310 and 330 degrees and between 240 and 270 degrees in Fig. 1 c). When the same analysis was performed on the data divided into four separate regional quadrants (Fig. 2) then a clearer directionality to road development became apparent, highlighting contrasting patterns of directional road development that were in general directed towards the centre of the Amazon, along the arc of deforestation, although a dominant direction was not apparent in the south-west quadrant (Fig. 2).

Figure 2 Radial plots of directional anisotropy of the Amazon divided by four quadrants. (a) North west (NW), (b) north east (NE), (c) south east (SE), (d) south west (SW). The direction in which development is moving is much more pronounced when four regions are considered as opposed to the Amazon as a whole (Fig. 1).

Model selection

For the travelling wave models, predictive performance wasmarginallybetter when fitted separately to the data divided into regional quadrants, explaining on average 25 % (95% CI = 12–35%) of the variance, than when the same model was fitted to all of the data together (explaining on average 22%; 95% CI = 12–29%). The parameter estimates (Appendix 1, Fig. S2, see supplementary material at Journals.cambridge.org/ENC) for the wave models corroborate the findings of the spatial anisotropy analysis; with the dominant direction of road development tending to be perpendicular to the arc of deforestation. The results also indicate that the average speed of wave movement was c. 54 km yr⁻¹, with the south-east having the highest speed of 114 km yr⁻¹ and the north-west having the lowest speed of23 km yr⁻¹. Given these speeds, the wave of road development would take on average 55 years (min = 26, max = 130 years) to traverse the study region (the region is approximately 3000 km wide along its widest dimension). This estimate is close to, and overlaps the bounds of, the wave speed calculated for the entire Amazon of 45 km yr⁻¹ (95% CI = 15–76).

The model performance metrics clearly indicate the exponential model with no neighbourhood effects (Eq. 3) to be the worst performing model in terms of its likelihood and ability to predict the evaluation data (Fig. 3). In contrast, the logistic model with no neighbourhood effects (Eq. 2) only appears to perform marginally worse than the models with neighbourhood effects, as indicated by a slightly higher DIC, with a mean DIC of 413 compared to mean DICs between 374 and 383, and lower training log likelihoods (TL), with a median TL of –0.8 compared to median TLs of between –0.74 and –0.75. However the correlation coefficient (CC) and the coefficient of determination (CD) indicate a very similar level of predictive accuracy to the models with neighbourhood effects.

Figure 3 (a) Goodness of fit measures for logistic (logi), exponential (exp) and fourneighbourhood effects models (Eqs 4–7). Mean parameter values and 95% confidence intervals are displayed for all goodness of fit measures except DIC, for which mean DIC (dark circles) and 10 DICs (grey circles) from each of the 10 subset parameter estimations are displayed. DIC = deviance information criterion, CC = coefficient of correlation, CD = coefficient of determination, TL = training likelihood, EL = evaluation likelihood. (b) Estimated parameters. Mean parameter values and 95% confidence intervals are displayed. D = magnitude of neighbourhood effect (units differ depending on formulation), K = maximum road density (km km⁻²), r = maximum road growth rate (km km⁻² yr⁻¹), τ = road density threshold difference (between neighbours) at which neighbourhood effects become apparent (km km⁻²), theta = estimated variance in the observations about the model predictions.

The best performing neighbourhood effects models predicted the log change of road density in the evaluation datasets with a correlation coefficient of c. 0.45 and a coefficient of determination of c. 0.2 (Fig. 3), the latter implying that the model explains c. 20% of the variance in the data. There are only minor quantitative differences in the predictive performance of the different neighbourhood effects models (Fig. 3). The maximum local rate of change in road density, r, tends to be higher for the logistic model with no neighbourhood effects than the neighbourhood effects models (Fig. 3; see Discussion for an explanation).

Comparing the predictions of road densities using the best wave and neighbour models (marginally the best in the case of the neighbourhood effects models) highlights contrasting abilities to predict log road densities and log changes in road densities. Comparison of predicted and observed log road densities for these models indicated a much better predictive performance by the neighbourhood model (Fig. 4). The NEm1 model predictions were closely correlated with the observed road densities (Spearman's correlation coefficient, r = 0.98, 95% CI = 0.97–0.98, p = < 0.001, t = 79.0, df = 283), while the wave model had a lower, but still significant, correlation (r = 0.54, 95% CI = 0.46–0.62, p < 0.01, t = 10.9, df = 283). Comparison of the two models in their ability to predict the log change in road density over the three-year period again revealed contrasting performances, while the neighbourhood model still had closer correlation between observed and predicted change in road density (r = 0.65, 95% CI = 0.58–0.71, p < 0.01, t = 14.5, df = 283, compared to r = 0.32, 95% CI = 0.21–0.42, p < 0.01, t = 5.7, df = 283). In the case of predicting change in road density, the wave model tended to over-predict the change in road density at low road densities, whereas the neighbourhood model was less biased (Fig. 4 b), although it still tended to over-predict log changes at low road densities.

Figure 4 (a) Observed road density in 2007 and road density in 2007 predicted by NEm1 and Amazon-wide wave models. (b) Observed versus predicted log_e road density in 2007 from wave (blue circles) and NEm1 (green circles) based on average median estimates for each location. Correlation lines for each model are displayed (solid lines, wave = blue, NEm1 = green). Correlations for upper and lower 95%confidence intervals are also displayed (dashed lines). A 1:1 line is shown for reference (red solid line). NEm1 has better predictions of absolute log_e road density in 2007. (c) Assessment of model predictive accuracy based on observed and predicted log_e density change between 2004 and 2007 (same colour scheme as in b).

Future projections

In terms of temporal scale, projections using both the neighbourhood effects andwave models were similar, with road density taking 65 years (based on 90% of grid cells reaching a road density greater than 0.47 km/km²) and c. 55 years, respectively, to reach the maximum predicted density across the entire Amazon (Fig. 5). Given that the assumption of no barriers to road development is unrealistic, wealsoincorporated barriers to road development (rivers and protected areas) on a 100-km grid in a simple way (where any grid cell with an area of more than 75% covered by barriers was considered ‘protected’, i.e. no roads would develop; Fig. 5). The results indicated very little influence of the barriers on the rate of road spread; only a few large areas were unaffected by road development, although this is most likely due to the coarse spatial resolution we adopted. We then repeated the analysis on a 10-km grid by simulating the model as a partial differential equation (for further details see Fig. 5). In this case, large areas of the Amazon remained relatively undeveloped and the rate of advance of road density was notably slowed by the barriers to development.

Figure 5 (a) Future projections of road density modelled on a 100-km grid, based on NEm1 (Eq. 5). Estimates suggest that within c.60 years the whole Amazon will have a relatively homogenous road density of 0.5 km km⁻². (b) Projections of future road density, modelled on a 100-km grid, based on NEm1 incorporating barriers to dispersal (road development),rivers and protected areas. A grid cell is considered a barrier when > 75% of its area is covered by a barrier (river or protected area). (c) Roads modelled on 10-km grid; the rate of spread is slowed and more complicated patterns of road development are evident when compared to projections made without real world dispersal barriers and at coarser resolutions. Simulation results in (a) and (b) are the mean estimates from running simulations for all 10000 combinations of parameter values for each of the 10 subset model fitting runs. The simulations in (c) were made by converting the model to a partial differential equation combining dispersal and logistic population growth, where the dispersal rate was determined by the NEm1 model formulation. It was solved using a fully explicit finite-difference method (Smith 1986) with time step of 0.01 years for one draw of parameters from one of the Markov Chains, simply to provide a representative simulation that would illustrate the effects of barriers at finer resolution.

DISCUSSION

Road network development forms an important part of regional development in Brazil, with the Brazilian government encouraging infrastructure development including the construction of new roads and paving of existing roads (Carvalho et al. 2002). The exact nature of this development will have influences on the spatial spread of roads. For example, the early development in the Amazon from the 1960/1970s and the newer Initiative for the Integration of the Regional Infrastructure of South America (IIRSA 2013) development plans focus on building long connection roads between areas, while the development plans of the 1990/2000s focused on improving existing road structures. Both approaches are intended to integrate the Amazon through colonization and development of roads, agriculture and industry, simultaneously boosting the economy and increasing living standards. However, they will have differing influences on the rate of road spread, with plans focusing on long roads increasing the rate at which the network traverses the Amazon. There are many factors that influence the spatio-temporal development of roads, besides direct government investment (Geist & Lambin 2002; Koorey 2009). However, despite the complex socioeconomic drivers behind road development our results imply that the emergent process at larger spatial scales can be approximated by a simple logistic growth and dispersal process. These findings support those of Ahmed et al. (2013), who, using empirical observations, determined that road development displayed a logistic growth pattern through time which appears to indicate that roads follow a ‘boom and bust’ dynamic.

We found little difference in model predictive performance when using different neighbourhood effects models so we avoided interpreting why we obtained subtle differences in predictive performance for those different dispersal functional forms. However all neighbourhood effects models supported unidirectional neighbourhood influences, from the region of high road density to the region of low density. This is not surprising, and similar neighbourhood effects are seen with other anthropogenic phenomena at large spatial scales, for example economic development (Sachs & Warner 1997) and land-use land-cover change (Verburg et al. 2004). The development of roads away from the arc-of-deforestation, towards the centre of the Amazon reflects the economic activities of the area, where initial roads grant access to extractive industry and colonizers, who expand the network with unofficial roads to increase access and transport products (Fearnside 2008; Perz et al. 2008). Over time, as more timber and agricultural land resources are exhausted, roads are built to access forest further from the arc-of-deforestation. Moving away from the arc-of-deforestation also helps integrate remote Amazon regions; an aim of initiatives such as Avanca Brazil (Fearnside 2008). This invites the question whether there are any other similar spatio-temporal patterns in the region that may drive, correlate with, or be driven by road development. Certainly, boom and bust dynamics have been observed in economic development of commodities within the Brazilian Amazon (Godfrey 1990; Macedo & Anderson 1993; Clough et al. 2009) and it is known that roads are a key driver of development (Chomitz & Gray 1996; Laurance et al. 2009). However, further analyses are required to assess quantitatively the degree to which waves of economic development and road development are related, in addition to looking for other qualitative and quantitative associations with road development at large spatial scales. We also infer that roads will tend to develop faster if the surrounding neighbours have high road densities. This is consistent with what is understood from patterns of economic development; for example Sachs & Warner (1997) found that countries which have more prosperous neighbours tend to do better economically, and in deforestation, whereareas with higher deforestation occur in close proximity to already deforested locations (Soler et al. 2007; Mann et al. 2010).

The estimated speeds of road development appear somewhat alarming, with both wave and neighbourhood effects models estimating complete road coverage across the Amazon within 60 years (year 2067), slightly faster than the 75 years estimated by Ahmed et al. (2013) but within the bounds of their estimates (58–103 years). Estimates suggest that by 2050, 40% of the Amazon forest will be lost (Soares-Filho et al. 2006) and that c. 44% of the forest will be under some form of protection (Verissimo et al. 2011). Such estimates imply that the majority of unprotected forest will be deforested by 2050, consistent with our projections. However our projections should still interpreted only as what would happen if the inferred pattern of road development between 2004 and 2007 were to continue for the next 60 years. While the estimated rate of 55 km yr⁻¹ is likely to continue in the short term, long-term rates remain uncertain. The fact that the model was only calibrated to changes in road density over three years, together with the relatively low predictive accuracy of the model (explaining about 20% of the variation in the data), mean that these projections should be treated as ‘what if’ scenarios (in this case, answering: ‘What if road density was to change in time and space exactly as projected by the model?’). Clearly multiple factors will influence the rate of road development to prevent this from happening, such as barriers to road development (as illustrated in Fig. 5), but also factors such as human population growth, economic development and trade policies, not considered in the models presented here.

The maximum local rate of change in road density, r, tends to be higher for the logistic model with no neighbourhood effects than the neighbourhood effects models. This is clearly due to higher overall local changes in road density having to be inferred for the logistic model to explain the observed changes in the absence of neighbourhood effects models. Similar maximum road densities are inferred for the logistic model with no neighbourhood effects and all of the neighbourhood effects models (Fig. 3). The superior predictive performance of the neighbourhood effects model compared to the wave model is largely because its projection incorporates the actual road density in 2004 (it only predicts the change over that time window, compared to the wave model, which simply predicts the log road density in 2007).

Our best fitted models explain around 20% of the variation in the data, clearly leaving room for improving predictions by incorporating more mechanisms. This would be a natural area for future work and could be conducted by extending the models we developed here within our parameter inference methodology (our source code can be downloaded from http://research.microsoft.com/en-us/downloads/b0cf61db-3c9d-4154-b5c1--5e5f72655185/default.aspx). One direction would be to infer the effects of different barrier types on the rate of road development. We performed a preliminary investigation into this by inferring parameters to a neighbourhood effects model, in which we reduced the maximum road density in a given cell (parameter K) in proportion to the area of that region that would act as a barrier to road development (in this case, rivers and protected areas). Protected areas are seen as being effective ways to protect habitats from anthropogenic destruction (Naughton-Treves et al. 2005; Nepstad et al. 2006; Joppa et al. 2008), however inclusion of protected areas in our preliminary investigation led to no detectable improvement in the model fit or the inferred parameters. Subsequent detailed investigation of the empirical data indicated many areas where roads developed in high density along river edges and extended into protected areas. Barreto et al. (2006) showed that road development occurred at a rate three times slower within protected areas compared to the surrounding area, but that it did still occur, and other studies have highlighted that protected areas offer far from perfect protection (Andam et al. 2008; Maiorano et al. 2008; Nagendra 2008). These, together with our preliminary investigation imply that that the existence of protected areas and barriers may have no discernible effect on the rate of secondary road development at large spatial scales. These results imply that barriers do get transgressed, and this is more likely to happen when the surrounding road density becomes high. This highlights an obvious area for future work. Incorporating and formalizing the relationships between road development and protected areas, indigenous reserves and other barriers, could enable more realistic future projections. Another avenue would be investigate and incorporate road ‘attractors’. For example areas of high agricultural suitability are likely to have greater road development susceptibility. These options could be explored by proposing process-based functional relationships between the rate of change of road density and these predictor variables and estimating the empirical support for them in the road density data as we have done in this study.