Supersampling and Network Reconstruction of Urban Mobility

Understanding human mobility is of vital importance for urban planning, epidemiology, and many other fields that draw policies from the activities of humans in space. Despite the recent availability of large-scale data sets of GPS traces or mobile phone records capturing human mobility, typically only a subsample of the population of interest is represented, giving a possibly incomplete picture of the entire system under study. Methods to reliably extract mobility information from such reduced data and to assess their sampling biases are lacking. To that end, we analyzed a data set of millions of taxi movements in New York City. We first show that, once they are appropriately transformed, mobility patterns are highly stable over long time scales. Based on this observation, we develop a supersampling methodology to reliably extrapolate mobility records from a reduced sample based on an entropy maximization procedure, and we propose a number of network-based metrics to assess the accuracy of the predicted vehicle flows. Our approach provides a well founded way to exploit temporal patterns to save effort in recording mobility data, and opens the possibility to scale up data from limited records when information on the full system is required.


Introduction
The increased pervasiveness of information and communication technologies is enabling the tracking of human mobility at an unprecedented scale. Massive call detail records from mobile phone activities [1,2] and the use of global positioning systems (GPS) in large vehicle fleets [3] for instance, are generating extraordinary quantities of positional and movement data available for researchers who aim to understand human activity in space. Other data sources, such as observations of banknote circulation [4,5], online location-based social networks [6,7], radio frequency identification traces [8][9][10], or even virtual movements of avatars in online games [11] have also been used as proxies for human movements. These studies have provided valuable insights into several aspects of human mobility, uncovering distinct features of human travel behavior such as scaling laws [4,12] or predictability of trajectories [13] among others. Besides empirical studies, the surge of available data on human mobility has also evoked interest in developing new theoretical models of mobility at several scales. Such models have deep implications for various subjects ranging from epidemiology to urbanism [14][15][16][17], with special importance in city planning and policy action [18].
Despite these first success stories, the theoretical development of tools and techniques for handling massive data sets of human mobility and for assessing their possible biases is still a road full of obstacles. Existing models based on gravity [19], radiation [20], intervening opportunities [21], etc. present a first step towards an accurate proxy for mobility at medium and large range scales, but they have been proven to be not always satisfactory to describe short scale movement such as intra-city displacements. The size of the data analyzed, the multiple scales involved, the highly skewed statistical nature of human activities [22] and the lack of strict control on the reliability of the data are just some of the multiple challenges this exciting new era poses.
Although they are often extensive, one of the main limitations of data sets used in empirically driven urban-scale mobility research is the limited coverage of the entire population under study. For instance, cell phone data records are typically obtained from a single operator. Similarly, data from taxis, or from other vehicle fleets are typically obtained from a single company, which usually represents only a small fraction of the actual number of vehicles circulating in a city [3,23]. In some scenarios, fully grasping a certain mobility-related phenomenon may require modelling the entire population of interest. For instance, it was shown that the fraction of taxi trips that can be shared in the city of New York is an increasing (albeit not simple) function of the number of daily taxi trips [24]. Hence, if a certain data set covers only a fraction of the daily taxi trips performed in a city, the taxi sharing potential cannot be fully unveiled.
The above discussion motivates the need of extrapolating urban mobility data starting from a subset of the population of interest. Although a number of urban mobility studies have applied such methods [25,26], a definition and assessment of a statistically rigorous extrapolation methodology is so far lacking. Even the sub-problem of assessing the quality of urban movement models is to date open, since the skewness of the underlying statistical distributions [10] makes a set of consistent, quantitative indicators hard to develop. In this paper, we fill these gaps by introducing a rigorous methodology to tackle the problem of obtaining an accurate picture of a mobility process when only a limited observation of such a process is available, both in time and volume. We first propose a simple rescaling rule which allows to quantify the strong temporal regularity of urban mobility patterns, even at very fine scales such as trips between particular intersections. Exploiting this regularity, we use a maximum entropy approach combining empirical data to model the occurrence of the core of frequent trips with an exponential gravity model [27][28][29] accounting for the variation observed in the least-frequent trips. We apply our method to accurately reconstruct the data set of all taxi trips performed in the city of New York in the year 2011 using small fractions sub-sampled from only a month of recorded data. By analysing the temporal patterns and the topological properties of the yearly mobility of taxis represented as a multi-edge network, we can finally assess the statistical accuracy of the proposed supersampling methodology using a number of both information-theoretical and networkbased performance metrics.
The remainder of the paper is structured as follows: We first present the study of both the temporal and topological patterns observed in the data, which then allows us to construct a maximum entropy method that exploits these features to solve the supersampling problem. Finally, we systematically test our reconstruction model on a very large data set and conclude by discussing some insights about the structure of urban mobility that the present study draws.

Results
Typically, mobility data is formalized by so-called Origin-Destination (OD) matrices, which are particular examples of weighted, or multi-edge networks [30]. OD matrices represent the number of observed tripst ij between the L = N 2 pairs of N locations or nodes i, j over a given observation period τ. A location can be defined based on a spatial partitioning of the urban area, on points of interest [31], or on road intersections [24]-as it is the case in the NY taxi data set at hand (see methods). Given this network representation, one can compute the total incomingŝ in j ¼ P it ij and outgoing strengthŝ out i ¼ P jt ij of a node i. Throughout this paper, we define active nodes as the subset of nodes which are either origin or destination of at least one trip in the set of all recorded tripsT ¼ P ijt ij and similarly active edges as the pairs of locations with at least one trip (t ij ! 1) recorded between them. The notationx shall refer to the observed value of the random variable x as derived from the data set, hxi to its expected value over independent realizations of a given model, while x denotes the matrix or network average of the variablex across the full empirical OD matrix (for example, average graph-degree k). Finally, the symbol hxi τ is used to express averages over time of variable x using bins of temporal length τ.

Stability of temporal urban-mobility patterns
While the built structure of cities evolves slowly in time, many dynamic, behavioral processes that take place within a city unfold relatively fast, and in principle could be strongly variable across time. However, human activity in cities exhibits highly regular patterns when observed over well defined periods of time, such as circadian or weekly rhythms. Intra-urban mobility is a good example for such activities: With longer time spans or larger samples of gathering movement data in cities, the picture of the underlying mobility network will clear up continually, but stable patterns should already emerge with relatively few data points as we can see in Fig 1A. To systematically test this hypothesis, we make use of a fleet of taxis acting as probes, sampling from the total traffic of all vehicles in a city. The total number of recorded trips, or sampling sizeT ðtÞ, depends on the total observation time τ and the number of probes, i.e., the size of the sub-population that is being monitored.
The evolution of the sampling size of trips as a function of the observation period, Fig 1B, can be extremely well approximated by a linear relation (R 2 > 0.999), indicating that the total number of trips generated daily in the city can be described as a random variable strongly concentrated around its mean value hTi τ = 1 day ' 403,000 ± 61,000 (confidence bounds reported as standard deviation).
On a yearly scale, the distribution of trips per day T τ = 1 day is not statistically compatible to a Gaussian distribution mainly due to seasonal effects, see Fig 1A. The effects of summer and winter holidays are apparent, and of Hurricane Irene that hit New York city towards the end of August, but if we disregard such outliers, corresponding to around 5% of the data that lies further than two standard deviations away from the mean, the quantile-quantile plot shows acceptable agreement with a normal distribution, see Fig 1C. To observe whether this strong regularity is also present in the finer structure of mobility, we must focus on each of the N nodes and L intersection pairs. Yet, we must find a suitable scaling to the data: The accumulated observed strength (both incoming and outgoing)ŝ i ðtÞ of each node and the weightt ij ðtÞ of each intersection pair will increase as more and more data is gathered, but if we normalize their (in or out) strength and weight by the total number of observed trips in the period τ, a strong regularity is recovered as we show in the following. The quantitiesp out;in intersection pair compared to the overall network, are extremely stable as shown in Table 1. We have split our data set into n τ equal time intervals (on daily, weekly and four-week bases) and computed the relative dispersion of the values accumulated over the entire data setpðt ¼ t max ¼ 1 yearÞ around the measured values hpi t , where τ max is the time at the end of the full observation period and the averages are performed over all the time slices of length τ. The graph-average of ε is very close to zero and highly doi:10.1371/journal.pone.0134508.g001 Table 1. Variability of node and node-pair statistics (incoming p in s and outgoing p out s relative strength, and relative number of trips between intersections p ij ) using different temporal granularity of the data set averaged over the full network (nodes and pairs of nodes respectively) compared to final yearly values. Time units with a total number of trips at least two standard deviations apart from the adjusted yearly mean have not been considered in the average to account for seasonal variations (< 5% of the data in the worst case). For the pairs of intersections ij, only pairs with at least one non-zero appearance on the time slicing have been considered for the average. The fraction of data with absolute relative error larger than two standard deviations is also reported as Outliers.
Time windows n τp in s : concentrated around this value for all the time windows considered (with a standard deviation of 13% in the worst case, decaying as sampling time is increased). Fig 2 shows the correlation between the relative error and the relative importance of nodes and links. The fact that coupled with second order seasonality effects induces an uneven distribution of errors: An overestimation of some values in the collection {hp ij i} will forcefully induce an underestimation in some other values of the collection. Despite this issue, we can clearly see that the vast majority of the mass of relative errors is concentrated around zero (see points in background for Fig 2).
To some extent, we would expect the node strength to be stable over time, since the number of trips received and generated at each location depends on parameters such as population density, number of points of interest present in a given location, etc. [32], whose evolution is given by much slower dynamics than the mobility process studied herein. But additionally, time stability is also observed at the trip level between intersections, yet in this case the analysis displays a higher variability around the mean-see Fig 2A and Table 1. This higher variability can be explained by the different sampling processes: While the percentage of active nodes becomes extremely stable already when just a very small number of days is considered (Fig 1B), this is not the case for the total number of active edges, because sampling of edge-specific attributes requires to grow as L ' N 2 to achieve a comparable level of accuracy.

The anatomy of urban flows
The results on temporal stability indicate that any model aiming to reproduce human mobility at urban scales should consistently exhibit regularities as reported above. Having seen that the mobility patterns are stable across time, an understanding of the main topological aspects of the aggregated static picture is further needed in order to be able to select the main features our methodology should aim to reproduce. The most relevant topological aspect of the mobility network is the highly skewed concentration of taxi pick-ups and drop-offs across the city, which gives rise to heavy-tailed node-strength distributions, Fig 3A (other general metrics are reported in the methods section). Therefore, to test whether this relevant property alone already captures the essential features of the mobility network, we must consider a null model which randomizes the considered network keeping constant the strength of each node -the multi-edge configuration model (see methods section or [30] for more details). When comparing the null model with the data, we observe significant deviations showing the importance of certain places or nodes in the network. In other words, even the strong heterogeneity of the distribution of strengths cannot account for the skewness of the weight distribution: Additional factors add many more trips between some connections than there should be under random conditions. The empirical link weight distribution, blue line in Fig 3B, is more skewed than under a random allocation of trips (Configuration model), dashed line, and the connections at the node level show a clear assortative correlation instead of a flat profile, see Fig 3D. This occurs despite the fact that the average number of trips between the most busy locations can be characterized by the relation t / s out i s in j with a reasonable accuracy, see Fig 3C, coinciding with the configuration model. These insights indicate that the distribution of node strengths across the city has a strong influence on the topology of the network, and needs to be taken into account when modelled, but needs extra ingredients to fully account for the observed pattern of connections.
A flexible model to reproduce human mobility A general maximum entropy based theory for model generation (see [27,29,33] for extended discussion and references) allows us to efficiently exploit both the observed temporal stability features and the heterogeneous topological properties of the network to solve the supersampling problem at hand. It starts with the assumption that each intersection pair in the network is allocated a constant fraction p ij of the total sampling (Eq (1)). Under this condition, and assuming that the mobility process is driven by some general constraints, such as population density or budget, it can be proved that for any desired level of sampling T d the statistics of trips for each pair of nodes can be well described by a set of L independent Poisson processes with mean ht ij i = T d hp ij i, with Following the theory, it would seem clear that from knowing the real values of the collection {hp ij i}, supersampling a mobility data set would be a trivial operation of generating L independent Poisson processes using the provided proportionality rule. Therefore, the problem now reduces to inferring the collection of values {p ij } from an available data set. We shall assume that only one snapshot of the aggregated mobility network is available to this end (thus assuming no temporal information is available on the trip data) as is usually the case in mobility studies. The maximum likelihood estimation of such values corresponds to There is, however, a practical issue in this formula related with the normalization condition for the random variables {p ij } and the presence of empty intersection pairs in the available observed data. For such intersections, using the formulas above, we have thatp ij ¼t ij ¼ 0, whereas their real hp ij i value is unknown but fulfils hp ij i 2 ½0;p min $T À1 . Since by definition both collections {p ij } and fp ij g need to be normalized, and denoting the set of active edges as from which we see that P ijj2E hp ij i P ijj2Epij ¼ 1. Hence, in general we cannot consider the empirically observed probabilitiesp ij as a good proxy for the real values of p ij unless the number of empty intersection pairs is very reduced. Given that the percentage of active edges (pairs of nodes for whicht ij > 0) is a very slowly increasing function of the sampling, see Fig 1B, inferring directly the set of probabilities {p ij } empirically would take an enormous data set-note that even with over a year of data only roughly 40% of edges are covered.
For the reasons given above, a simple proportionality rule using Eq (4) is not a good supersampling strategy, specially for skewed and sparse data sets.

Supersampling urban trips
Based on the previous discussion, we now present the methodology for supersampling an urban mobility data set that consists in inferring the collection of {hp ij i} values from a set of aggregated empirical tripst ij n o . We should do so bearing in mind that despite the maximum likelihood formula in Eq (4) cannot be directly used for the empty intersection pairs in the data, it does perform well for non-empty intersections (see Fig 2). The maximum entropy based framework naturally allows such a procedure, since it can combine any constraint driven model with the rich information encoded in the trip sample. We propose a method to predict trips based on the theory mentioned earlier: Taking the L intersection pairs (being them active edges in the data set or not), we split them into two parts, the subgroup of trusted trips defined as Q ¼ ij jt ij > t min n o and its complementary part Q C .
The value t min is a threshold modelling a minimal statistical accuracy that depends on the amount of data available, and which may be set to 1 in practical applications. We keep the proportionality rulep ij /t ij for the trusted trips, while for the remaining trips we apply a doubly constrained exponential gravity model -other maximum entropy models [29] could also perform well as long as they preserve the outgoing and incoming strength. In other words, we generate a collection of {hp ij i} values, 8 > <

> :
The values γ and {x i , y j } are the 2N + 1 Lagrange multipliers satisfying the following equationŝ whereĈ ¼ P ijt ij c ij is the total euclidean distance of the observed trips (c ij stands for the distance between intersections i and j). Note that, by construction, the values are properly normalized, i.e., ij2Q C x i y j e Àgc ij ¼ 1. The model presented earlier needs to deal with the issue of inactive nodes that do not appear in the original data due to poor sampling, i.e., nodes for whichŝ ¼ 0 either incoming or outgoing for some observation period τ. This issue has a minor impact in our case due to the previously observed rapid coverage of the number of active nodes (the number of inactive nodes is negligible after accumulating very few days of data, see Fig 1B). In any case, it can be solved easily: given that the geographic positions of the nodes are available, we could always artificially assign a certain relative strength to the nodes not present in the data using complementary call detail records [34], census data or points of interests (POI) data, or assign them some values according to a chosen distribution depending on the data at hand. For simplicity, in our case we have chosen to keep only the nodes present in the original data.

Assessing the quality of the supersampling methodology
To test the supersampling methodology, we have proceeded to select a timespan of our data set corresponding to an observation period of τ = 1 month (February 2011) from which we further randomly sub-sample different fractions f used as training sets to compute {hp ij i} applying Eqs (6) and (7). We then reconstruct the OD using the proportionality rule {ht ij i = T d hp ij i} for both the complete and reduced data set, T d = T(τ 0 = 1 year) and T d = T(τ = 1 month). Finally, we compare the model predictions with the set of empirically observed trips in these periods.
In order to do so, we need to introduce metrics to quantify the resemblance between model predictions and actual recorded values. Commonly used indicators are the Sorensen-Dice common part of commuters (CPC) value [35] or the linear fit of ht ij i model vst ij [20]. However, the skewness of the observed trip distribution represents a challenge to these indicators: The variability of low-valued trips induces notable instabilities on both and being single numbers, they are only able to provide a limited picture on the precision of a given model to reproduce empirical results.
To overcome these issues, we propose a slight modification to both the Sorensen CPC index and the coefficient of determination R 2 from the fit ht ij i model vst ij (see methods) and we additionally propose the introduction of a number of network metrics to precisely assess the quality of the models used in human mobility at the topological, finer, scale: The unweighted degrees of the nodes hkiðŝÞ, the weighted neighbor strength correlation hs w nn iðŝÞ, existing trip distribution P(t) and number of existing trips as a function of the origin destination strength product tðss 0 Þ. Table 2. Parameters for the validation of the methodology. See details on each indicator in the main text and in the methods section. The number of trusted trips fed to the model relative to the entire number of generated trips f Q ¼ P ijjij2Qt^ij =T d ðtÞ is reported in column 3. The Supersampled models with different fractions f are only generated using subsamples of the training set (1 month observation period). Empirical stands for the model generated using the empirical probabilitiesp^i j (Eq (1)) of the full data set and Configuration stands for the multi-edge configuration model applied to the full data set.  The results for the supersampling method are summarized in Table 2 and a specific example for f = 0.1 (reconstruction using only 10% of the original data of the monthly data set compared to yearly data) is shown in Fig 4 for  We observe an accurate reconstruction of the mobility network for a wide range of values of f, which shows the validity of our proposed supersampling methodology. At the global scale, even at extreme levels of subsampling, our model is successful at reconstructing the original dataset. Also at the topological scale, despite the heterogeneities in the underlying distributions, the methodology generates very accurate predictions. The predictions for the least frequently visited nodes display higher relative errors due to the presence of inactive nodes in the training dataset (1.6% of total nodes for f = 0.1).
Upon close inspection, our inferred values {hp ij i} slightly over-estimate low-valued weights and underestimate large-valued weights and strengths as was expected from the analysis in temporal stability (see Fig 2), yet the errors are greatly mitigated as we can see in Fig 3B. See Figs 3B (green line) and 4E (green dots) where we can observe a gap around t * 100 and ht þ ij i $ 100, respectively, which corresponds to the separation between the trusted empirical data (separated points in the background belonging to the group of trusted trips Q) and the reconstructed trips (clustered cloud of points). The minor seasonal fluctuations detected first in our temporal analysis together with these over-and under-estimations explain the minor limitations of the model to reproduce perfectly the entire yearly data set.
The second order effects induced by the seasonality of recorded data can also be seen in the performance of our methodology under extreme levels of subsampling (using around 1% of the sample monthly data to feed the model or less). In these circumstances, the model is still able to produce a good prediction of the empirical data, yet it reproduces better the accumulated yearly mobility rather than the monthly one since the inherent seasonal variations of traffic between certain intersection pairs are smoothed by the aggregation procedure.
Furthermore, in the event that enough historical data were available, we could achieve even better results by computing the collection {hp ij i τ } with an appropriate τ period (depending on the granularity of the data) and approximating hp ij i ' hp ij i τ for the group of trusted trips (Eq (7)). Such a procedure, which may be extended to overcome the minor limitations imposed by the seasonality of the data and other improvements related with the presence of non-active nodes could be derived to perfect the method.

Discussion
The stationarity of the temporal statistics of trips between different locations, together with a suitable scaling for the data observed over different timespans, has allowed us to develop a general model that is highly effective at reconstructing a general mobility scenario from very limited aggregated data, without the need of having fine grain temporal information, and provides insights about the structure of real taxi trips. The success of our reconstruction method, even using very small amounts of data, points out the composite structure of the network of urban mobility: Taxi displacements are characterized by a small core of very frequent trips coupled with trips generated at random but conditioned by the structural constraints of the city such as population distribution and mobility costs.
Ultimately, the results presented in this paper could be used to answer questions that are of fundamental importance in the field of human mobility modelling, such as: i) Can an accurate picture of urban mobility patterns be obtained from an incomplete sample of the population?, and ii) are existing metrics sufficient to assess the quality of model predictions?
We pointed out the importance of data sampling and of the correct assessment of mobility models, and introduced network-based tools to evaluate such models. The implications of these findings are two-fold: On the one hand, the stationarity of the temporal patterns could be exploited to save space and effort in recording mobility data. On the other hand, our method opens the possibility of efficiently scale up data from reduced fleet of vehicles in cases where a full knowledge of the system is needed.
Our study provides a first step in showing that incomplete samples can indeed be scaled up adequately with the appropriate models, and that network metrics are required to comprehensively assess mobility model predictions.

Materials and Methods
Data set OD matrices are typically inferred from either census/survey data or alternative means such as social media data or Call Detail Records (CDR) [36]. They must be constructed in terms of trips, i.e. well defined trajectories between a starting and ending point. For this reason, we use the taxi data for its completeness, since all trips are recorded, and we consider it a good proxy for general urban mobility, with the caveat that commuting patterns and areas with little taxi demand are covered less. For the present paper we have used the full set of taxi trips (with customers) with starting and ending points within Manhattan obtained from the New York Taxi and Limousine Commission for the year 2011 via Freedom of Information Law request [24]. We have aggregated the data at the node level taking into account the grid of roads up to secondary level and adding trips to the nearest node (intersection) in a radius of 200m. We have decided to keep the self-loops present in the data for simplicity (albeit their fraction is completely negligible). In the analysis, all trips including both week-ends and week-days are considered, since the pattern for weekly trips shows a continuous increase in the number of trips peaking on Friday and followed by a sudden drop on Sundays.
For the subsampling of the reduced data set used as basis to reconstruct the networks, we have used random subsampling. The parameters of the obtained mobility networks for two different observation periods τ are reported in Table 3.
The analyzed data shows that most of the taxis share similar performance. This, coupled with the fact [25] that individual taxi mobility traces are in large part statistically indistinguishable from the overall population, justifies that their individual traces (corresponding to sets of trips performed by different customers which can be considered as independent events) can be safely aggregated for the analysis.

Null model: The Multi-Edge Configuration model
Cities usually display a high level of variability across its different locations in terms of activity, i.e., city center concentrate busy locations while outskirts usually display less traffic/retails areas and others. At the level of networks, this translates in important topological heterogeneities which need to be accounted for using a suitable null model. The Multi-Edge configuration model [30] is a maximum entropy model that is used to generate maximally random network instances of graphs with a prescribed strength sequence {s in , s out } (number of trips entering and leaving each node). The expected number of trips between two nodes i and j with respective strengths s out i and s in j reads ht Conf ij i ¼ŝ out iŝ in j =T and the assortativity profile is flat (nodes are uncorrelated at the level of strengths). Throughout this paper, we use this null model as benchmark preserving the strength sequence obtained from aggregating the complete yearly observation period.
Indicators for the quality of the reconstruction Distance based measures.
• Sorensen-Dice common part of commuters index: This indicator was proposed in [35] and based on the original formulation is defined as We propose an alternative version formulated in terms of averages which reads, The different versions of this indicator have values in the range [0, 1], where CPC = 1 indicates total coincidence between data and model and CPC = 0 total disagreement. However, for sparse data sets with a skewed distribution oft ij n o values, Eq (8) may return values excessively lower than 1, even for models very close to reality. To exemplify this fact, Table 4 shows a comparison of the performance of the two indicators for the models presented in Table 2. For the Empirical model, we can see that the second version of the indicator recovers values close to 1 as would be expected. Furthermore, both indicators converge to very similar values as sampling is increased (the yearly data set contains roughly 12 times more trips than the monthly one). • Linear correlation ht ij i model vs t^i j : This method is widely used [20]. We report the coefficient of determination R 2 cond in all tables, based on the comparison between real data and conditional values of the model on the existing edges (since we are using a biased statistic based only on the observed trips in the original data, not the entire set of intersection pairs). With Poisson distributed variables with mean ht ij i, such conditioned average reads, which converges to the average value for highly used trips. Hence, the coefficient of determination R 2 cond assuming an identity relation ht þ ij i ' t data is explicitly, Network measures. To better grasp the quality of the models, we propose to compare also some of its multi-edge network related quantities: • Degree: The unweighted degree of the nodes is the sum of their incoming/outgoing active edges (edges for which t ij > 0), k x = ∑ x Θ(t ij ) being E = ∑ ij Θ(t ij ) the total number of active edges and x referring to the outgoing i or incoming direction j. • Distribution of weights on existing edges: This commonly used measure is computed as P(t) = ∑ ij δ t,t ij /E and indicates the collection of weight values present in the network, where δ x,y corresponds to the Kroenecker delta.
• Graph average existing weight of trips as a function of product of incoming and outgoing degree: To quantify the deviation from a completely randomized configuration model, we compute also this metric, which is the average weight of existing trips as a function of the product of out(in) strengths of their origin (destination) nodes: tðss 0 Þ ¼ P ijjj¼s 0 ;i¼st ij =n ss 0 , where n ss 0 is the cardinality of the sum. For the configuration case, this magnitude is equivalent to ht þ i / ss 0 1Àe Àss 0 =T . Information values. We also assess the quality of our models using their Log-Likelihood values assuming a set of independent Poisson random variables with known means ht ij i for each intersection pair, Incompatible Loglikelihood values are not reported in tables (such as cases where ht ij i ¼ 0 < t ij or jt ij À ht ij i j) 0).

Simulations
All the simulations and solving of saddle point equations as well as the analysis of the multiedge networks have been performed using the freely available, open source package Origin-Destination Multi-Edge analysis (ODME) [37].