Indirect evidence of network effects in a system of cities

Abstract

We describe a simple spatial model of urban growth for systems of cities at the macroscopic scale, which combines direct interaction between cities and an indirect effect of physical network flows as population growth drivers. The model is parametrized on population data for the French system of cities between 1831 and 1999, in which strong non-stationarity in correlation patterns suggest to apply the model on local time windows. The corresponding calibration of the model using genetic algorithms provides the evolution of interaction processes and network effects in time. Furthermore, the fit improvement when adding network module appears effective when controlling for additional parameters, what confirms the ability of the model to unveil network effects in the system of cities.

Keywords

Urban systems urban growth spatial interactions network effects empirical Akaike Information Criterion

Introduction

Cities are paradoxically unsustainable and source of negative externalities, but also the best chance to reach sustainability and resilience to climate change (Glaeser, 2011). The dynamics of urban systems at a macroscopic scale, and more precisely drivers of urban growth, are crucial processes that need to be understood to meet these potentialities. A better knowledge of how cities differentiate, interact and grow is thus a relevant topic both from a theoretical perspective and for policy applications. Indeed, Pumain et al. (2009) suggest that cities are incubators of social change, their fate being closely linked to the one of societies.

Various disciplines have studied models of urban growth with different objectives and taking diverse aspects into account. For example, economics are still reluctant to include spatial interactions in models (Krugman, 1998) but are extremely detailed on market processes, even for models in economic geography. On the other hand, geography focuses more on territorial specificities and interactions in space but will have more difficulties to produce general conclusions. The example of these two disciplines and their misunderstandings (Marchionni, 2004) shows how it is difficult to make bridges, how it needs exceptional efforts to translate from one to the other (as P Hall did for Von Thunen work (Taylor, 2016)), and therefore how it is far from evident to grasp the complexity of urban systems in an integrated way.

The simplest model to explain urban growth, the Gibrat model, assumes random growth rates for cities. It has been shown by Gabaix (1999) that this model asymptotically produce the expected rank-size law (Zipf’s law) for system of cities, which is considered as one of the most regular stylized facts, at least in its generalized scaling law formulation (Nitsch, 2005). Explaining urban scaling laws is closely related to the understanding of urban growth, as Bettencourt et al. (2008) suggest that these reflect underlying universal processes and that all cities are scaled version of each other. This approach however does not reflect the complex relation between economic agents for which Storper and Scott (2009) advocate. Using a bottom-up reconstruction of urban areas with dynamical microscopic population data, Rozenfeld et al. (2008) show indeed that positive deviations to the rank-size law systematically exist, and that these must be an effect of spatial interactions between urban areas.

Complexity approaches are good candidates to integrate these interactions into models. Andersson et al. (2006) introduce for example a model of urban economy as a growing complex network of relations. The Evolutive Urban Theory, introduced by Pumain (1997), focuses on cities as co-evolving entities and produces explanations for growth at the scale of the system of cities. Pumain et al. (2006) show that scaling laws could be a consequence of functional differentiations and of the diffusion of innovation between cities. The positioning regarding universality of laws is more moderate than theories of scaling (West, 2017), as Pumain (2012b) highlights that ergodicity can difficultly be assumed in the frame of complex territorial systems. One crucial feature of this paradigm is the importance of interactions between agents, generally the cities, to produce the emergent patterns at the scale of the system. Pumain and Sanders (2013) investigate the advantages of agent-based models in comparison to more classical equation systems, and this methodological aspect is in accordance with the theoretical positioning, as it allows to take into account the heterogeneity of possible interactions, the geographical particularities, and to naturally include emergence between levels and render multi-scale patterns.

In this paper we aim at exploring further the assumption, central to Pumain’s Evolutive Urban Theory, that spatial interactions between cities are significant drivers of their growth. More precisely, we consider both abstract interactions and flow interactions mediated through the physical networks, mainly the transportation networks. We extend existing models accordingly. Our contribution is twofold: (i) we show that very basic interaction models based on population only can be fitted to empirical data and that adjusted parameter values can be directly interpreted; and (ii) we introduce a novel methodology to quantify overfitting in models of simulation, as an extension of Information Criteria for statistical models, in which application to our calibrated models confirms that fit improvement is not only due to additional parameters, but that the extended model effectively capture more information on system processes. This unveils network effects in an indirect way. We first review modeling approaches to urban growth based on spatial interactions.

Urban growth and spatial interactions

First of all, we must precise that we consider only models at the macroscopic scale, ruling out the numerous and rich approaches at the mesoscopic scale, that include for example cellular automata models, models of urban morphogenesis or land-use change models. Some include explicitly spatial interactions: Rybski et al. (2013) show that a gravity-based aggregation of population in a lattice produces realistic metropolitan patterns. Such kind of models operates however at an other level of detail than the one we consider. We also naturally rule out economic models that do not explicitly include spatial interactions.

Several models of urban growth at the macro scale have insisted on the role of space and spatial interactions. Frasco et al. (2014) introduce a spatialized model of social network growth that provides a good agreement with observed data both for city size hierarchy and the relation between population and density. Bretagnolle et al. (2000) propose a spatial extension of the Gibrat model. The gravity-based interaction model that Sanders (1992) uses to apply concepts of synergetics to cities is also close to this idea of interdependent urban growth, contained physically in the phenomenon of migration between cities.

A more refined extension with economic cycles and innovation waves is developed by Favaro and Pumain (2011), yielding a version of the core of Simpop models (Pumain, 2012a) in terms of dynamical systems. This family of models has started with a toy-model based on economic interactions between cities as agents that yield hierarchical patterns at the scale of the system (Sanders et al., 1997). The Simpop2 model, still based on spatial interaction for commercial exchanges, including successive innovation waves, unveils structural differences between the European and the US urban systems (Bretagnolle and Pumain, 2010). The SimpopLocal model (Pumain and Reuillon, 2017a) simulates the emergence of initial settlement patterns. The Marius model (Cottineau, 2014) couples population and economic growth with interactions between cities, allowing to accurately reproduce real trajectories on the former Soviet Union after calibration with processes multi-modeling.

Urban growth and transportation networks

Under similar assumptions of previously reviewed models, the inclusion of transportation networks has been rarely pursued, contrary to the mesoscopic scale at which relations between networks and territories have been widely studied by LUTI models for example (Chang, 2006). Transportation network growth models (Xie and Levinson, 2009), prolific in economics and physics, are designed to reproduce network topology and do not focus mainly on the influence of networks on cities. Bigotte et al. (2010) study an optimization model for network design combining the effects of urban hierarchy and of transportation network hierarchy. Baptiste (1999) models the dynamical interplay between network links capacity and city growth on a subset of the French city system. The SimpopNet model (Schmitt, 2014) goes a step further in modeling the co-evolution between cities and transportation networks, as it allows new network links to be created in time. These examples show the difficulty of coupling these two aspects of urban systems in models of growth, and we will for this reason take into account network effects in a simplified way as detailed further.

The rest of this paper is organized as follows: our model is introduced and formally described in the next section; we then describe the results obtained through exploration and calibration of the model on data for French cities, in particular the unveiling of network effects significantly influencing growth processes thanks to a novel methodology introduced. We finally discuss implications of these results.

Model description

Rationale

Some confusion may arise when surveying stochastic and deterministic models of urban growth. To what extent is a proposed model “complex” and is the simulation of stochasticity necessary? Concerning the Gibrat model and most of its extensions, independence assumptions and linearity produce a totally predictable behavior, which is thus not complex in the sense of exhibiting emergence, in the sense of weak emergence (Bedau, 2002). In particular, the full distribution of random growth models can be analytically determined at any time (Gabaix, 1999), and in the case of studying only the first moment, a simple recurrence relation avoids to proceed to any Monte-Carlo simulation. Under these assumptions, it is natural to work with a deterministic model, as it is done for example for the Marius model (Cottineau, 2014). We will work with that hypothesis, capturing complexity through non-linearity.

We work on simple territorial systems assumed as regional city systems, in which cities are basic entities. The time scale corresponds to the characteristic scale associated to this spatial scale, i.e. around one or two centuries. Spatial interactions are captured through gravity–type interactions, this simple formulation having the advantage of being simple and of capturing the first law of Tobler, namely that interaction strength fades with distance. Other approaches introduced recently perform similarly at this scale (Masucci et al., 2013).

Model description

We consider a deterministic extension of the Gibrat model, what is equivalent to consider only expectancies in time. Let $\vec{P} (t) = (P_{i} (t))_{1 \leq i \leq n}$ be the population of cities in time. Under Gibrat independence assumptions, we have $Cov [P_{i} (t), P_{j} (t)] = 0$ . A linear extended version would write $\vec{P} (t + 1) = R \cdot \vec{P} (t)$ where R is an independent random matrix of growth rates (identity in the original case). This yields directly thanks to the independence assumption that $E [\vec{P} (t + 1)] = E [R] \cdot E [\vec{P}] (t)$ . We generalize this linear relation to a non-linear relation that allows to be more consistent with model interpretation and more flexible. Denoting $\vec{μ} (t) = E [\vec{P} (t)]$ , we take

\vec{μ} (t + 1) = Δ t \cdot f (\vec{μ} (t))

(1)

Note that in that case, stochastic and deterministic versions are not equivalent anymore, precisely because of the non-linearity, but we stick to a simple deterministic version for the sake of simplicity. The specification of the interdependent growth rate is thus given by

f (\vec{μ}) = r_{0} \cdot Id \cdot \vec{μ} + G (\vec{μ}) \cdot \vec{1} + \vec{N} (\vec{μ})

(2)

where

\vec{1}

is the column vector full of ones, and

G = G_{ij} = w_{G} \cdot \frac{V_{ij}}{< V_{ij} >}

such that the interaction potential V_ij follows a gravity–type expression given by, with d_ij distance between i and j (Euclidian or network distance)

V_{ij} = (\frac{μ_{i} μ_{j}}{(\sum_{k} μ_{k})^{2}})^{γ_{G}} \cdot exp (- d_{ij} / d_{G})

(3)

The network effect term $\vec{N}$ is given by $N_{i} = w_{N} \cdot \frac{W_{i}}{< W_{i} >}$ where the network flow potential W_i reads

W_{i} = \sum_{k < l} (\frac{μ_{k} μ_{l}}{(\sum_{j} μ_{j})^{2}})^{γ_{N}} \cdot exp (- d_{kl, i} / d_{N})

(4)

where

d_{kl, i}

is the distance of city i to the shortest path between k, l computed in the geographical space, which can be through a transportation network or in an impedance field of the Euclidian network. All seven model parameters are detailed below.

The first component corresponds to the pure Gibrat model, that we obtain by setting the weights $w_{G} = w_{N} = 0$ . The second component captures direct interdependencies between cities, under the form of a separable gravity potential such as the one used in Sanders (1992). The third term is aimed at capturing network effects by expressing a feedback of network flow between cities k, l on the city i. Intuitively, a demographic and economic flow physically transiting through a city or in its surroundings is expected to influence its development (e.g. through intermediate stops), this effect being of course dependent on the transportation mode since a high speed line with few stops will skip most of the traversed territories. Note that we don’t use exactly gravity flows in the network term, since there is no decay of interactions generating flows with distance, but a decay of the effect of the flow as a distance to the network: it is equivalent to assuming a long-range use of the network on average in time, and is this way complementary to the first gravity term.

Model parameter space

We give in Table 1 the description of model parameters, detailing the associated processes and parameter ranges. Both direct interaction and second-order network flows effect have the same structure, namely separability between the effect of distance and the influence of population, an exponential decay parameter, and a hierarchy parameter expressing the inequality of contribution depending on relative sizes of cities: the highest the exponent, the more contribution of smaller cities will be negligible regarding larger cities.

Table 1.

Model parameters’ summary.^a

Parameter	Notation	Process	Interpretation	Range	Unit
Growth rate	r ₀	Endogenous growth	Growth rate	$[0, 1]$	1
Gravity weight	w_G	Direct interaction	Max average rate	$[0, 1]$	1
Gravity gamma	γ_G	Direct interaction	Level of hierarchy	$[0, + \infty]$	1
Gravity decay	d_G	Direct interaction	Interaction range	$[0, + \infty]$	km
Feedback weight	w_N	Flows effect	Max average rate	$[0, 1]$	1
Feedback gamma	γ_N	Flows effect	Level of hierarchy	$[0, + \infty]$	1
Feedback decay	d_N	Flows effect	Network effect range	$[0, + \infty]$	km

We give the parameters’ names, notations, associated processes, possible interpretations, typical variation ranges, and units.

We propose to interpret the distance decay parameter in the following way. Let us fix an arbitrary fraction α and typical spatial ranges for a local urban system d_L and for a long range urban system d_R, consider a city i and two neighbors $j, j'$ with same population $μ_{j} = μ_{j}'$ , at distances d_L and d_R of i, respectively. If we want to answer the question to what distance difference is equivalent to an attenuation of α of the interaction potential with i, we obtain $d_{L} - d_{R} = - d_{G} \cdot \ln α$ . Therefore, d_G is exactly the proportionality coefficient answering this intuitive request.

Finally, we will consider only positive weights w_G and w_N, to follow empirical observations as detailed below. Numerical values for the weights will be given normalized by number of cities implied in the process, i.e. $w'_{G} = w_{G} / n$ and $w'_{N} = w_{N} / (n (n - 1) / 2)$ .

Data

Our model is assumed as hybrid as it relies on semi-parametrization on real data. It could be possible to study it as a full toy-model, initial configuration and physical environment being constructed as synthetic data. We however aim at unveiling stylized facts on real data rather than on model behavior in itself, and setup therefore the model from the data we now describe.

Population data

We work with the Pumain-INED historical database for French cities (Pumain and Riandey, 1986), which give populations of Aires Urbaines (INSEE definition) at time intervals of five years, from 1831 to 1999 (31 observations in time). The latest version of the database integrates Urban Areas, allowing to follow them on a long time-period, according to the long time ontology for cities given by Bretagnolle (2009), that constructs a functional definition of cities as entities with boundaries evolving in time. We work on the 50 bigger cities in 1999. We furthermore isolate periods of similar length excluding wars, obtaining nine periods of 20 years¹ on which adjustment of the model will be done, taking into account non-stationarity in time.

Physical flows

The model takes into account stylized physical flows, in order to be rather independent of network effective shape. These are therefore assumed to follow the geographical shortest path according to terrain slope. This assumption avoids to obtain geographical absurdities such as cities with a difficult access having an overestimated growth rate. Using a 1 km resolution Digital Elevation Model, we construct an impedance field, following Collischonn and Pilar (2000), which is given by

Z = (1 + \frac{α}{α_{0}})^{n_{0}}

where Z is the impedance of links the 1 km grid network in which each cell is connected to its eight neighbors. α is the terrain slope computed with elevation difference between the two cells. We take fixed parameter values

α_{0} = 3

(corresponding to approximatively the real world value of a 5% slope) and

n_{0} = 3

which yield more realistic paths than smaller values.²

Indicators of model performance

We work on an explanatory rather than an exploratory model. Therefore, indicators to evaluate model outputs are not directly linked to intrinsic properties of trajectories or obtained final states, but rather to a distance to the phenomenon we want to explain, i.e. the data. Given real population $p_{i} (t)$ (historical realizations of $P_{i} (t)$ ) and simulated expected populations $μ_{i} (t)$ obtained with $\vec{μ} (t_{0}) = \vec{p} (t_{0})$ on a period of length T, we can evaluate two complementary aspects of model performance:

• Overall model performance, given by logarithm of the mean-square error in space and time

ɛ_{G} = \ln (\frac{1}{T} \sum_{t} \frac{1}{n} \sum_{i} (p_{i} (t) - μ_{i} (t))^{2})

• Average local model performance, given by the mean-square error on logarithms

ɛ_{L} = \frac{1}{T} \sum_{t} \frac{1}{n} \sum_{i} (\ln p_{i} (t) - \ln μ_{i} (t))^{2}

Both are actually complementary, as using only $ɛ_{G}$ as it is generally done will focus only on larger cities and give poor results on medium-sized and small cities (for France only Paris will have reasonable fit as it strongly dominates other urban areas and cities). $ɛ_{L}$ allows therefore to take into account model performance in all cities simulated by the model.

Results

Stylized facts

Basic stylized facts can be extracted from such a database, as it has already been widely explored in the literature (Guérin-Pace and Pumain, 1990). We retrieve better fits of log-normal distributions of growth rates at all dates compared to normal fits, and also the fact that growth rates are mainly positive, on the cities we consider and when removing wars.

An interesting feature to look at in relation with our considerations on spatial interactions is correlations between time-series, and more particularly their variation as a function of distance. We consider 50 years overlapping time-windows to have enough temporal observations, finishing respectively in 1881, 1906, 1931, 1962, 1999, and estimate on each, for each couple of cities (i, j), the correlation between log-returns ${\hat{ρ}}_{ij} = ρ [Δ X_{i}, Δ X_{j}]$ with a classical Pearson estimator, where

Δ X_{i} = X_{i} (t) - X_{i} (t - 1)

and

X_{i} (t) = \ln (\frac{P_{i} (t)}{P_{i} (t_{0})})

This method, used mainly in econophysics (Mantegna and Stanley, 1999), reveals dynamical interactions without being biased by sizes.

We show in Figure 1 the smoothed correlations curves as a function of distance, for each time period. First of all, the strong differences between each time period confirm the non-stationarity of growth rates over the whole time period, and justify the use of local fit in time for the model. We can also interpret these patterns in terms of historical events for the system of city and the transportation network. The dynamic of the system begins with a flat correlation in 1881, around 0.2, that could be spurious due to simultaneous similar growth for all cities. It then stays flat but goes to zero, witnessing strong differentiations in growth patterns between 1881 and 1906. After 1931, the effect of the distance is clear with decreasing curves, starting between 0.4 and 0.5. We postulate that this evolution must be partly linked to the evolution of the transportation network: considering railway network for example (Thévenin et al., 2013), the initial overall development may have fostered long range interactions flattening thus the correlation curves, whereas its maturation over time has conducted to the return of more classical interactions decreasing quickly with distance.

Figure 1.

Time–series correlations as a function of distance. Solid line corresponds to smoothed correlations, computed between each pairs of normalized log-returns of population time-series, on successive periods given by curve color.

Model exploration

Implementation

Data preprocessing, result processing and models profiling are implemented in R. For performance reasons and an easier integration into the OpenMole software for model exploration (Reuillon et al., 2013), a scala version was also developed. The question of a trade-off between implementation performance and interoperability is a typical issue in this kind of model, as a fully blind exploration and calibration can be misleading for further research directions or thematic interpretations. A NetLogo implementation, allowing interactive exploration and dynamical visualization, was also developed for this reason. Source code for models, cleaned raw data, simulation data, and results used in this paper are available on the open repository of the project at https://github.com/JusteRaimbault/InteractionGibrat.

We show in Figure 2 an example of model output. Cities color gives city-level fit error and their size the population. Outliers can therefore easily be spotted (as Saint-Nazaire having the worst fit in the example shown) and possible regional effects identified. We illustrate in pink an example of geographical shortest path, from Rouen to Marseille, which reasonably corresponds to the actual current shortest path by highway. Top right plot shows trajectory in time for a given city, whereas the bottom right plot gives overall fit quality in time, by plotting simulated data against real data. The closest the curve is to the diagonal, the better the fit.

Figure 2.

Example of output of the model. The graphical interface allows to explore interactively on which cities’ changes operate after a parameter change, what is necessary to interpret raw calibration results. The map gives adjustment errors by city (color) and their population (size). We illustrate in pink the geographical shortest path between Rouen and Marseille. The plot in top right panel gives in time the trajectory of a selected city, comparing simulated population with real population. The bottom right plot gives for each date the simulated data against real data: the closest the curve is to the diagonal, the better the fit.

Behavior patterns

First model explorations, by simply sweeping fixed grids of the parameter space, already suggest the presence of network effects, in the sense that physical flow effectively has an influence on growth rates. We show in Figure 3 a configuration of such a case. At fixed gravity parameters and growth rate, we study variations of the parameters w_N, d_N, and γ_N and the corresponding response of $ɛ_{G}$ and $ɛ_{L}$ . At fixed values of γ_N, we observe similar behaviors of the indicators when w_N and d_N change. The existence of a minimum for both as a function of d_N, that becomes stronger when w_N increases, shows that introducing the network feedback terms improves local and global fits compared to the Gravity Model alone, i.e. that the associated process have potential explanatory power for growth patterns.

Figure 3.

Evidence of network effects revealed by model exploration. Left plot gives $ɛ_{G}$ as a function of d_N for varying $r_{0} / w_{N}$ , at fixed gravity effect and $γ_{N} = 3$ . Right plot is similar for $ɛ_{L}$

Calibrating the Gravity Model

We now use the model to indirectly extract information on processes in time. Indeed, under assumption of non-stationarity, the temporal evolution of locally adjusted parameters shows the evolution of the corresponding aspect of the processes.

In a first experiment, we set w_N = 0 and calibrate the model with four parameters on the nine successive 20-year time windows. The optimization problem associated to model calibration does not present features allowing an easy solving (closed-form of a likelihood function, convexity or sparsity of the optimization problem), we must rely on alternative techniques to solve it. Brute force grid search is rapidly limited by the dimensionality curse. Classical methods (Batty and Mackie, 1972) such as gradient descent fail because of the rather complicated shape of the optimization landscape. Calibration using genetic algorithms (GA) is an efficient solution to find approximate solutions in a reasonable time. OpenMole embeds a collection of such meta-heuristics for different purposes: Schmitt et al. (2015) demonstrate the capabilities of these methods to calibrate models of simulation. In our case, it furthermore allows to do a bi-objective calibration on $(ɛ_{G}, ɛ_{L})$ .

We use the standard steady state GA provided by OpenMole (Pumain and Reuillon, 2017b), distributed on 25 islands, with a population of 200 individuals and 100 generations. We show in Figure 4 the calibration results on successive periods, by plotting final populations in the objective space. As expected, Pareto fronts that correspond to compromises between the two opposite objectives are the rule. It means that the model cannot be accurate both globally and locally, and an intermediate solution has to be found. This may due to the fact that interaction range changes with city size (i.e. that terms in the potential are no longer separable), that we keep as a possible model development. The shapes of the Pareto front are revealing the chaotic optimization landscape, as for some periods such as 1921–1936 or 1962–1982 fronts are not regular and sparse. The change in shapes also translates different dynamical regimes across the periods: for 1881–1901, the quasi-vertical shape followed by an isolated front at high $ɛ_{G}$ values reveals a quasi-binary model behavior in the optimal regimes, in the sense that improving $ɛ_{L}$ under the limit is only possible through a qualitative jump at a high price for $ɛ_{G}$ . The values taken by d_G for periods 1892–1911 and 1921–1936 show that larger cities have longer interaction ranges, as higher values of d_G give better values of $ɛ_{G}$ .

Figure 4.

Calibrating the Gravity Model. Pareto-front on successive periods. Each panel gives the Pareto front between the two objectives $(ɛ_{G}, ɛ_{L})$ , the dates for the period being given at the top of each. Color level gives the value of distance decay parameter d_G.

We show in Figure 5 the values of fitted parameters in time, averaged over the Pareto front and for best single-objective solutions. First, the two peak patterns for r₀ correspond roughly to the patterns observed in average growth rates. The evolution of w_G has a similar shape but lagged by 20 years: it can be interpreted as a repercussion of endogenous growth on interaction patterns in the following years, which is consistent with an interpretation of the interaction process in terms of migrations. The values of d_G, with an increase until 1900 followed by a progressive decrease, is consistent with the behavior of empirical correlations commented above: the first 50 years windows in Figure 1 have a higher interaction range corresponding to flat correlation curves. Finally, the level of hierarchy γ_G has regularly decreased, corresponding to an attenuation of the power of large cities that can be understood in terms of the progressive decentralization in France that has been fostered by the administration.

Figure 5.

Calibrated parameters for Gravity Model only. Each plot gives fitted values in time t for each parameter: in order from left to right and top to bottom, growth rate r₀, gravity weight w_G, gravity decay d_G, gravity gamma γ_G. Red and Green curves correspond to best points for $ɛ_{G}$ (respectively $ɛ_{L}$ ), whereas the blue curves give the average value over the Pareto front with standard deviation.

Unveiling network effects

We now turn to the calibration of the full model on successive periods, in order to interpret parameters linked to network flows and gain insight into network effects. The full calibration is done in a similar way with seven parameters being free. We plot in Figure 6 the fitted values in time for some of these parameters.

Figure 6.

Calibrated parameters for the full model. We plot values of $r_{0}, w_{G} / r_{0}, d_{N}$ and γ_N in time t, for single-objective optimal points (Red and Green curves) and averaged over the Pareto front (Blue).

The behavior of growth rate and of the gravity weight relative to growth rate, that is similar to the Gravity Model only, confirms that network effects are indeed at the second order and that endogenous growth and direct interactions are main drivers. Network effects are however not negligible, as they improve the fit as shown before in model exploration, capturing therein second-order processes. The evolution of d_N, corresponding to the range on which network influences the territories it goes through, shows a minimum in 1921–1936 to stabilize again later, but at a value lower than past values. This could correspond to the “tunnel effect”, when high-speed transportation does not stop in territories it goes through. Indeed, the evolution of railway has witnessed a high decrease for local lines at a date similar to the minimum, and later the emergence of specific high speed lines, explaining this lower final value. Hierarchy of flows has slightly decreased as for gravity, but is extremely high. This means that only flows between larger cities have a significant effect. This way, the model gives indirect information on the processes linked to network effects.

Evaluating model performance

We focus in this last experiment on quantifying the “performance” of the model, taking into account not only its predictive capabilities, but also its structure. More precisely, we want to tackle the issue of overfitting, which has been for long recognized in the field of machine learning for example (Dietterich, 1995), but for which there is a lack of methods for models of simulation. We need to introduce a tool to confirm that the improvement in model fit is not only artificially due to additional parameters.

The Akaike Information Criterion (AIC) provides, for statistical models for which a likelihood function is available, the gain in information between two models (Akaike, 1998), correcting fit improvement for number of parameters. Similar methods include the Bayesian Information Criterion (BIC), which relies on slightly different assumptions. Biernacki et al. (2000) propose an integrated likelihood as a generalization of these criteria for unsupervised classification. Pohle et al. (2017) show that in the case of selecting the number of states in Hidden Markov Models, real cases induces too much pitfalls for standard methods to work robustly, and suggest pragmatic selection based on their results and expert judgment. In our case, the problem is that it is not even possible to define the criteria.

The method we propose is based on the intuitive idea of approaching models of simulation by statistical models and using the corresponding AIC under certain validity conditions. Bastani et al. (2017) use a similar trick of considering the models as black boxes and approaching them to gain insights, in their case to extract interpretable structure as decision trees.

Let (X, Y) be the data and observations. We consider computational models as functions $(X, α_{k}) \to M_{α_{k}}^{(k)} (X)$ mapping data values to a random variable. What is seen as data and parameters is somehow arbitrary but is separated in the formulation as corresponding dimensions will have different roles. We assume that the models have been fitted to data in the sense that an heuristic has been used to find an approximate optimal solution $α_{k}^{*} = {argmin}_{α_{k}} ‖ M_{α_{k}}^{(k)} (X) - Y ‖$ , and we write $ɛ_{k} = ‖ M_{α_{k}}^{(k)} (X) - Y ‖^{2}$ the corresponding mean-square error. For each optimized computational model, a statistical model $S^{(k)}$ with the same degree of freedom can be fitted on a set of realizations: $M_{α_{k}^{*}}^{(k)} (X) = S^{(k)} (X)$ , with an error $s_{k} = {‖ M_{α_{k}^{*}}^{(k)} (X) - S^{(k)} (X) ‖}^{2}$ .

We follow the intuition that if statistical models are good approximations of the simulation models compared to models discrepancy to reality, namely $s_{k} << ɛ_{k}$ , then the gain of information between the two should mostly capture the gain of information between simulation models.³ We define therefore an Empirical AIC measure between two simulation models by

I (M^{(1)}, M^{(2)}) = Δ AIC [S^{(1)}, S^{(2)}]

(5)

In practice we calibrate the gravity only model and the full model on the full time span, and choose two intermediate solutions giving $M^{(1)}$ at $r_{0} = 0.0133, d_{G} = 4.02 e 12, w_{G} = 1.28 e - 4, γ_{G} = 3.82$ with $ɛ_{G} = 31.2375, ɛ_{L} = 302.89$ and the full model M⁽²⁾ at $r_{0} = 0.0128, d_{G} = 8.43 e 14, w_{G} = 1.230 e - 4, γ_{G} = 3.81, w_{N} = 0.60, d_{N} = 7.47 e 14, γ_{N} = 1.15$ with $ɛ_{G} = 31.2366, ɛ_{L} = 302.93$ . It is not clear how the empirical method is sensitive to the type of statistical model used. We use therefore several models for robustness, each time with the corresponding number of parameters (4 for the first and 7 for the second model): a polynomial model of the form $a_{0} + \sum_{i > 0} a_{i} X^{i}$ , a mixture of logarithm and polynomial as $a_{0} + a_{1} \ln X + \sum_{i > 1} a_{i} X^{i}$ , and a generalized polynomial with real power coefficients that are optimized for model fit using a genetic algorithm $a_{0} + \sum_{i > 0} a_{i} X^{α_{i}}$ . We fit the statistical models using successive years as different realizations.

Results are shown in Table 2. We give the value of

s_{k} / ɛ_{k}

and the

Δ AIC

. We also provide the

Δ BIC

, that correspond to the Bayesian criterion differences on statistical models, to check the robustness regarding the information criterion used. We find a positive value for five criteria out of six, confirming that information gain is indeed positive. The gain decreases when statistical model fit improves, and only the BIC for the optimized model fails to show an improvement. The assumption of negligible errors is always verified as the rate stays around 1%. This approach is of course preliminary and further work should be done for a more systematic testing and a more robust justification of the method. It suggests however that fit improvement in the model of simulation is effective, and that the model reveals therefore network effects.

Table 2.

Empirical AIC results.

Statistical model	M⁽¹⁾ Relative fit	M⁽²⁾ Relative fit	ΔAIC	ΔBIC
Polynomial	0.01438	0.01415	19.59	3.65
Log-polynomial	0.01565	0.01435	125.37	109.43
Generalized polynomial	0.01415	0.01399	11.70	−4.23

Discussion

Multi-modeling

We obtain for most of periods and calibration points large values of γ_N in the calibration of the full model, what suggests that a model with only major flows could be a possible alternative. An additional parameter determining the cutoff value for flows would have to be introduced, and to what extent this variation improves model performance should be tested in future developments. More generally, this raises the issue of more systematic multi-modeling approaches and of benchmarking models. Cottineau et al. (2015) introduce a modular framework to compare concurrent hypothesis to explain urban growth. In the same spirit, an immediate extension of our work would be to test different functional specifications and benchmark more diverse network processes, such as the diffusion of innovation as done by Favaro and Pumain (2011) or the creation of social ties (Frasco et al., 2014).

Theoretical implications

Our results support the hypothesis that physical transportation networks are necessary to explain the morphogenesis of territorial systems, in the sense that some aspects are fully contained within networks and cannot be approximated by abstract proxies. We showed indeed on a relatively simple case that the integration of physical networks into some models effectively increases their explanative power even when controlling for overfitting. This can be understood as a direction to expand Pumain’s Evolutive Urban Theory (Pumain, 1997). This theory considers networks as carriers of interactions in systems of cities but does not put a particular emphasis on their physical aspect and the possible spatial patterns resulting from it such as bifurcations or network induced differentiations. The development of a sub-theory focusing on these aspects is an interesting direction suggested by our empirical and modeling results.

Urban system specificity

The model has not yet been tested on other urban systems and other temporalities, and further work should investigate which conclusions we obtained here are specific to the French urban system on this period, and which are more general to any system of cities. Applying the model to other system of cities also recalls the difficulty of defining urban systems. In our case, a strong bias should arise by considering France only, as Lille must be highly influenced by Brussels for example. The extent and scale of such models is always a delicate question. We rely here on the administrative coherence and the consistence of the database, but sensitivity to system definition and extent should also be further tested.

Towards co-evolutive models

Our focus on network effects remains quite limited since (i) we do not consider an effective infrastructure but abstract flows only, and (ii) we do not take into account the possible network evolution, due to technical progresses (Bretagnolle et al., 2000) and infrastructure growth in time. An interesting development would be first the application of our model with real network data, using effective distance matrices in time, computed, e.g. with the train network used by Thévenin et al. (2013).

Then, allowing the network to dynamically evolve in time, as a function of flows, would yield a model of co-evolution between cities and transportation networks for a system of cities. This co-evolution has been highlighted empirically by Bretagnolle (2009). This kind of model is not frequent, and Schmitt (2014) provides with the SimpopNet model, one of the few examples. It is shown by Raimbault (2017) that disciplinary compartmentalization may be at the origin of the relative absence of such type of models in the literature. Indeed, it would imply to include heterogenous processes such as economic rules to drive network growth that are quite far from the approach taken. It would however allow to investigate to what extent the refinement of network spatial structure and network dynamics can improve the explanation of urban system dynamics.

Conclusion

We have introduced a spatial model of growth for a system of cities at the macroscopic scale, including second-order network effects among endogenous growth and direct interaction growth drivers. The model is parametrized on real data for the French city system between 1831 and 1999. The calibration of the model in time provides interpretations for the evolution of processes of interaction within the system of cities. We furthermore show that the model effectively unveils network effects by controlling for overfitting. This work paves the way for more complicated models with dynamical networks that would capture the co-evolution between transportation network and territories.

Footnotes

Acknowledgements

Results obtained in this paper were computed on the vo.complex-system.eu virtual organization of the European Grid Infrastructure (). We thank the European Grid Infrastructure and its supporting National Grid Initiatives (France-Grilles in particular) for providing the technical support and infrastructure.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Notes

Juste Raimbault is a PhD candidate at Université Paris 7, UMR CNRS 8504 Géographie-cités and UMR-T -T IFSTTAR 9403 LVMT, under the supervision of Arnaud Banos and Florent Le Néchet. He holds an engineer degree from Ecole Polytechnique and Ecole des Ponts et Chaussães, and a Master in Complex Systems Science from Ecole Polytechnique. His research focuses on complexity approaches to the relations between networks and territories.

References

Akaike H (1998) Information theory and an extension of the maximum likelihood principle. In: Parzen E, Tanabe K and Kitagawa G (eds) (2012). Selected Papers of Hirotugu Akaike. New York: Springer-Verlag, pp. 199–213.

Andersson

Frenken

Hellervik

(2006) A complex network approach to urban growth. Environment and Planning A 38(10): 1941.

Baptiste H (1999) Interactions entre le système de transport et les systèmes de villes: Perspective historique pour une modélisation dynamique spatialisée. PhD Thesis, Centre d’études supérieures de l’aménagement (Tours).

Bastani O, Kim C and Bastani H (2017) Interpretability via model extraction. arXiv preprint arXiv:1706.09773.

Batty

Mackie

(1972) The calibration of gravity, entropy, and related models of spatial interaction. Environment and Planning A 4(2): 205–233.

Bedau

(2002) Downward causation and the autonomy of weak emergence. Principia: An International Journal of Epistemology 6(1): 5–50.

Bettencourt

Lobo

West

(2008) Why are large cities faster? Universal scaling and self-similarity in urban organization and dynamics. The European Physical Journal B: Condensed Matter and Complex Systems 63(3): 285–293.

Biernacki

Celeux

Govaert

(2000) Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7): 719–725.

Bigotte

Krass

Antunes

, et al.(2010) Integrated modeling of urban hierarchy and transportation network planning. Transportation Research Part A: Policy and Practice 44(7): 506–522.

10.

Bretagnolle A (2009) Villes et réseaux de transport: Des interactions dans la longue durée, France, Europe, États-Unis. Hdr, Université Panthéon-Sorbonne – Paris I.

11.

Bretagnolle A, Mathian H, Pumain D and Rozenblat C (2000) Long-term Dynamics of European Towns and Cities: Towards a Spatial Model of Urban Growth. Cybergeo: European Journal of Geography.

12.

Bretagnolle A and Pumain D (2010) Comparer deux types de systèmes de villes par la modélisation multi-agents. In: Weisbuch G and Zwirn H (eds) Paris : Vuibert, pp. 271–299.

13.

Chang

(2006) Models of the relationship between transport and land-use: A review. Transport Reviews 26(3): 325–350.

14.

Collischonn

Pilar

(2000) A direction dependent least-cost-path algorithm for roads and canals. International Journal of Geographical Information Science 14(4): 397–406.

15.

Cottineau C (2014) L’évolution des villes dans l’espace post-soviétique. Observation et modélisations. PhD Thesis, Université Paris 1 Panthéon-Sorbonne.

16.

Cottineau

Reuillon

Chapron

, et al.(2015) A modular modelling framework for hypotheses testing in the simulation of urbanisation. Systems 3(4): 348–377.

17.

Dietterich

(1995) Overfitting and undercomputing in machine learning. ACM Computing Surveys (CSUR) 27(3): 326–327.

18.

Favaro

Pumain

(2011) Gibrat revisited: An urban growth model incorporating spatial interaction and innovation cycles. Geographical Analysis 43(3): 261–286.

19.

Frasco

Sun

Rozenfeld

, et al.(2014) Spatially distributed social complex networks. Physical Review X 4(1): 011008.

20.

Gabaix X (1999) Zipf’s law for cities: an explanation. The Quarterly Journal of Economics 114(3): 739–767.

21.

Glaeser E (2011) Triumph of the City: How Our Greatest Invention makes us Richer, Smarter, Greener, Healthier, and Happier. Penguin.

22.

Guérin-Pace

Pumain

(1990) 150 ans de croissance urbaine. Economie et statistique 230(1): 5–16.

23.

Krugman

(1998) Space: The final frontier. The Journal of Economic Perspectives 12(2): 161–174.

24.

Mantegna RN and Stanley HE (1999) Introduction to Econophysics: Correlations and Complexity in Finance. Cambridge: Cambridge University Press.

25.

Marchionni

(2004) Geographical economics versus economic geography: Towards a clarification of the dispute. Environment and Planning A 36(10): 1737–1753.

26.

Masucci

Serras

Johansson

, et al.(2013) Gravity versus radiation models: On the importance of scale and heterogeneity in commuting flows. Physical Review E 88(2): 022812.

27.

Navakatikyan

(2007) A model for residence time in concurrent variable interval performance. Journal of the Experimental Analysis of Behavior 87(1): 121–141.

28.

Nitsch

(2005) Zipf zipped. Journal of Urban Economics 57(1): 86–100.

29.

Pohle J, Langrock R, van Beest F, et al. (2017) Selecting the number of states in hidden markov models – Pitfalls, practical challenges and pragmatic solutions. arXiv preprint arXiv:1701.08673.

30.

Pumain

(1997) Pour une théorie évolutive des villes. Espace géographique 26(2): 119–134.

31.

Pumain D (2012a) Multi-agent system modelling for urban systems: The series of Simpop models. In: Heppenstall AJ, Crooks AT, See LM and Batty M (eds) Agent-based Models of Geographical Systems. Netherlands: Springer, pp. 721–738.

32.

Pumain D (2012b) Urban systems dynamics, urban growth and scaling laws: The question of ergodicity. In: Portugali J, Meyer H, Stolk E and Tan E (eds) Complexity Theories of Cities Have Come of Age. Berlin: Springer, pp. 91–103.

33.

Pumain D, Paulus F and Vacchiani-Marcuzzo C (2009) Innovation cycles and urban dynamics. In: Lane D, Pumain D, van der Leeuw SE and West G (eds) Complexity Perspectives in Innovation and Social Change. Dordrecht: Springer, pp. 237–260.

34.

Pumain D, Paulus F, Vacchiani-Marcuzzo C and Lobo J (2006) An Evolutionary Theory for Interpreting Urban Scaling Laws. Cybergeo: European Journal of Geography.

35.

Pumain D and Reuillon R (2017a) The simpoplocal model. In: Pumain D and Reuillon R (eds) Urban Dynamics and Simulation Models. 22: Springer, pp. 21–35.

36.

Pumain D and Reuillon R (2017b) Urban Dynamics and Simulation Models. Springer International Publishing.

37.

Pumain

Riandey

(1986) Le fichier de l’ined. Espace, Populations, Sociétés 4(2): 269–277.

38.

Pumain

Sanders

(2013) Theoretical principles in interurban simulation models: A comparison. Environment and Planning A 45(9): 2243–2260.

39.

Raimbault J (2017) Models coupling urban growth and transportation network growth: An algorithmic systematic review approach. Plurimondi (17).

40.

Reuillon

Leclaire

Rey-Coyrehourcq

(2013) Openmole, a workflow engine specifically tailored for the distributed exploration of simulation models. Future Generation Computer Systems 29(8): 1981–1990.

41.

Rozenfeld

Rybski

Andrade

, et al.(2008) Laws of population growth. Proceedings of the National Academy of Sciences 105(48): 18702–18707.

42.

Rybski

Ros AGC Kropp

(2013) Distance-weighted city growth. Physical Review E 87(4): 042114.

43.

Sanders L (1992) Système de villes et synergétique. Economica.

44.

Sanders

Pumain

Mathian

, et al.(1997) Simpop: A multiagent system for the study of urbanism. Environment and Planning B 24: 287–306.

45.

Schmitt C (2014) Modélisation de la dynamique des systèmes de peuplement: de SimpopLocal à SimpopNet. PhD Thesis, Paris 1.

46.

Schmitt

Rey-Coyrehourcq

Reuillon

, et al.(2015) Half a billion simulations: Evolutionary algorithms and distributed computing for calibrating the simpoplocal geographical model. Environment and Planning B: Planning and Design 42(2): 300–315.

47.

Storper

Scott

(2009) Rethinking human capital, creativity and urban growth. Journal of Economic Geography 9(2): 147–167.

48.

Taylor PJ (2016) A polymath in city studies. In: Knowles RD and Rozenblat C (eds) Sir Peter Hall: Pioneer in Regional Planning, Transport and Urban Geography. Springer International Publishing, pp. 11–20.

49.

Thévenin

Schwartz

Sapet

(2013) Mapping the distortions in time and space: The French railway network 1830–1930. Historical Methods: A Journal of Quantitative and Interdisciplinary History 46(3): 134–143.

50.

West G (2017) Scale: The Universal Laws of Growth, Innovation, Sustainability, and the Pace of Lifein Organisms, Cities, Economies, and Companies. Penguin.

51.

Xie

Levinson

(2009) Modeling the growth of transportation networks: A comprehensive review. Networks and Spatial Economics 9(3): 291–307.