Supersampling and Network Reconstruction of Urban Mobility

Oleguer Sagarra; Michael Szell; Paolo Santi; Albert Díaz-Guilera; Carlo Ratti

doi:10.1371/journal.pone.0134508

Abstract

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.

Citation: Sagarra O, Szell M, Santi P, Díaz-Guilera A, Ratti C (2015) Supersampling and Network Reconstruction of Urban Mobility. PLoS ONE 10(8): e0134508. https://doi.org/10.1371/journal.pone.0134508

Editor: Matjaz Perc, University of Maribor, SLOVENIA

Received: April 20, 2015; Accepted: July 9, 2015; Published: August 14, 2015

Copyright: © 2015 Sagarra et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited

Data Availability: Data are freely available from the New York Taxi and Limousine Commission via Freedom of Information Law request. The dataset is composed by all the Taxi trips recorded in NY for the year 2011 by the NYTLM. Data requests may be sent to Jason Gonzalez, Records Access Officer: foil@tlc.nyc.gov.

Funding: The work of O.S. has been partially supported by the EU-LASAGNE Project, Contract No. 318132 (STREP), the Spanish MINECO Grant FIS2012-38266-C02-02 and by the Generalitat de Catalunya through the FI Program and 2014-SGR-00608. P.S. and C.R. thank the Enel Foundation, Audi Volkswagen, and all the members of the MIT Senseable City Lab Consortium for supporting this research. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: Audi Volkswagen partly funded this study. There are no patents, products in development or marketed products to declare. This does not alter the authors' adherence to all the PLOS ONE policies on sharing data and materials.

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–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–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–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 network-based 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 trips between the L = N² 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 and outgoing strength 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 trips and similarly active edges as the pairs of locations with at least one trip () recorded between them. The notation shall refer to the observed value of the random variable x as derived from the data set, ⟨x⟩ to its expected value over independent realizations of a given model, while denotes the matrix or network average of the variable across the full empirical OD matrix (for example, average graph-degree ). Finally, the symbol ⟨x⟩_τ 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 size , depends on the total observation time τ and the number of probes, i.e., the size of the sub-population that is being monitored.

Download:

Fig 1.

(a) Circadian clocks in the city: Time-series of the number of observed trips , active edges and active nodes in the city of NY per day over the year 2011. The influence of both seasonal fluctuations, major events and stable weekly patterns are clearly observed. (b) Aggregated fraction of active nodes and total trips as days of data are accumulated normalized by the total number of recorded nodes and trips at the end of the observation period, the evolution of the accumulated graph density in time E/L = E/N² (observed binary edges over total possible pairs of intersections) is also reported. Nodes are almost fully sampled within the first days analyzed, while edges are sampled sub-linearly in time. (c) Quantile—Quantile plot comparing the number of trips per day observed in the data set without outliers within two standard derivations of the mean (> 95% of the data) to the theoretical quantiles of a normal distribution and linear fit (dashed line) showing their proportionality (similar results not shown obtained for number of edges and nodes respectively).

https://doi.org/10.1371/journal.pone.0134508.g001

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² > 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 ⟨T⟩_{τ = 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) of each node and the weight 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 quantities and , which quantify the relative importance of a given node and intersection pair compared to the overall network, (1) 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 set around the measured values , (2) 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 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).

Download:

Table 1. Variability of node and node-pair statistics (incoming

and outgoing

relative strength, and relative number of trips between intersections

) 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.

https://doi.org/10.1371/journal.pone.0134508.t001

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 {⟨p_ij⟩} 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).

Download:

Fig 2. The effect of sampling on node and intersection pair temporal stability.

Correlation between measured values of intersection pair (a) and node statistics (b), (c) for different time slices and relative dispersion around the mean (Eq (2)) for the yearly aggregated data. Error bars represent standard deviations on the log-binned data. Raw data for the daily case is shown in the background. For visual clarity, panel a) only shows a random subsample of 1/1000 of the original points.

https://doi.org/10.1371/journal.pone.0134508.g002

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² 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 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.

Download:

Fig 3. Empirical network features of the taxi multi-edge mobility network.

(a) Distribution of incoming and outgoing strengths s. (b) Existing edge weight complementary cumulative distribution function compared with a configuration model [30] and the supersampled model with f = 0.1. (c) Weighted average neighbor out-strength as function of node out-strength (similar results for incoming strengths not shown). (d) Graph-average existing weight as a function of product of outgoing and incoming strengths of origin and destination for a single instance of the network. Points represent mean values and error bars represent standard deviations computed using log-binning and distributions also shown using log-binning and ensemble results averaged over 100 repetitions of the two models.

https://doi.org/10.1371/journal.pone.0134508.g003

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 ⟨t_ij⟩ = T_d⟨p_ij⟩, with (3)

Following the theory, it would seem clear that from knowing the real values of the collection {⟨p_ij⟩}, 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 (4) 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 that , whereas their real ⟨p_ij⟩ value is unknown but fulfils . Since by definition both collections {p_ij} and need to be normalized, and denoting the set of active edges as , we have, (5) from which we see that .

Hence, in general we cannot consider the empirically observed probabilities 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 which ) 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 {⟨p_ij⟩} values from a set of aggregated empirical trips . 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 and its complementary part 𝓠^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 rule 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 {⟨p_ij⟩} values, (6) The values γ and {x_i, y_j} are the 2N + 1 Lagrange multipliers satisfying the following equations (7) where 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., .

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 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 {⟨p_ij⟩} applying Eqs (6) and (7). We then reconstruct the OD using the proportionality rule {⟨t_ij⟩ = T_d⟨p_ij⟩} for both the complete and reduced data set, T_d = T(τ′ = 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 ⟨t_ij⟩^model vs [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² from the fit ⟨t_ij⟩^model vs (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 , the weighted neighbor strength correlation , existing trip distribution P(t) and number of existing trips as a function of the origin destination strength product .

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 the different indicators proposed. For comparison, results using both a configuration model and the empirical values using the information encoded in the complete dataset are also shown.

Download:

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 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 probabilities (Eq (1)) of the full data set and Configuration stands for the multi-edge configuration model applied to the full data set.

https://doi.org/10.1371/journal.pone.0134508.t002

Download:

Fig 4. Main network differences between real data accumulated over a year and Supersampled model from real one month data with f = 0.1 subsampling.

(a)-(d) Relative error ( with being a magnitude measured from the aggregated yearly network) between reconstructed network using supersampling and original data for outgoing degrees (a), strengths (b) and average neighbor strength (d) (similar results for incoming direction not displayed). The complementary cumulative distribution function of both edge lengths and trips lengths (c) is also shown. (e) Comparison between empirical values and model prediction over a single run. Configuration model expectation from a single run using the full year data set is also shown for comparison. All results averaged over 100 repetitions of the model, error bars represent standard deviations on the log-binned data and raw data is shown in the background. For visual clarity, panel e) only shows a random subsample of 1% of the raw data in background.

https://doi.org/10.1371/journal.pone.0134508.g004

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 {⟨p_ij⟩} 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 , respectively, which corresponds to the separation between the trusted empirical data (separated points in the background belonging to the group of trusted trips 𝓠) 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 {⟨p_ij⟩_τ} with an appropriate τ period (depending on the granularity of the data) and approximating ⟨p_ij⟩ ≃ ⟨p_ij⟩_τ 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.

Download:

Table 3. Empirical network parameters.

Symbol stands for graph averaged magnitudes. See main text and methods for a definition of each magnitude.

https://doi.org/10.1371/journal.pone.0134508.t003

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ⁱⁿ, s^out} (number of trips entering and leaving each node). The expected number of trips between two nodes i and j with respective strengths and reads 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 (8) We propose an alternative version formulated in terms of averages which reads, (9) 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 of 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 ⟨t_ij⟩^model vs : This method is widely used [20]. We report the coefficient of determination 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 ⟨t_ij⟩, such conditioned average reads, (10) which converges to the average value for highly used trips. Hence, the coefficient of determination assuming an identity relation is explicitly, (11)

Download:

Table 4. Values for the different versions of the common part of commuters index for the reconstructed models.

⟨CPC_sample⟩ and averaged trip values computed over 1000 repetition of each model, standard deviations lower than 10⁻³ for all cases. Note how differences between indicators disappear with increased sampling.

https://doi.org/10.1371/journal.pone.0134508.t004

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.
Average weighted neighbor strength: This metric is widely used in the literature. It indicates the level of correlations at the node level and is defined as , where y is the complementary of x (if x = i then y = j and vice versa).
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: , where n_ss′ is the cardinality of the sum. For the configuration case, this magnitude is equivalent to .

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 ⟨t_ij⟩ for each intersection pair, (12) Incompatible Loglikelihood values are not reported in tables (such as cases where or ).

Simulations

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

Acknowledgments

We thank P. Colomer-de-Simon for useful comments and suggestions.

Author Contributions

Conceived and designed the experiments: CR AD PS OS. Performed the experiments: OS. Analyzed the data: MS OS. Contributed reagents/materials/analysis tools: OS. Wrote the paper: AD PS MS OS.

References

1. González MC, Hidalgo CA, Barabási AL. Understanding individual human mobility patterns. Nature. 2008;453(7196):779–782. pmid:18528393
- View Article
- PubMed/NCBI
- Google Scholar
2. Blondel VD, Decuyper A, Krings G. A survey of results on mobile phone datasets analysis. arXiv preprint arXiv:1502.03406. 2015.
3. Bazzani A, Giorgini B, Rambaldi S, Gallotti R, Giovannini L. Statistical laws in urban mobility from microscopic GPS data in the area of Florence. J Stat Mech. 2010;2010:P05001
- View Article
- Google Scholar
4. Brockmann D, Hufnagel L, Geisel T. The scaling laws of human travel. Nature. 2006;439(7075):462–465. pmid:16437114
- View Article
- PubMed/NCBI
- Google Scholar
5. Thiemann C, Theis F, Grady D, Brune R, Dirk Brockmann D. The Structure of Borders in a Small World. PLoS ONE. 2010;5:e15422. pmid:21124970
- View Article
- PubMed/NCBI
- Google Scholar
6. Scellato S, Noulas A, Lambiotte R, Mascolo C. Socio-spatial Properties of Online Location-based Social Networks. Proc Int AAAI Conf Weblogs Soc Media. 2011;11.
- View Article
- Google Scholar
7. Scellato S, Musolesi M, Mascolo C, Latora V, Campbell A. NextPlace: A Spatio-Temporal Prediction Framework for Pervasive Systems. Pervasive’11 Proceedings of the 9th international conference on Pervasive computing. 2011;152–169.
8. Barthélemy M. Spatial Networks. Phys Rep. 2010;499:1–101.
- View Article
- Google Scholar
9. Cattuto C, Van den Broeck W, Barrat A, Colizza V, Pinton JF, Vespignani A. Dynamics of person-to-person interactions from distributed RFID sensor networks. PloS ONE. 2010;5(7):e11596. pmid:20657651
- View Article
- PubMed/NCBI
- Google Scholar
10. Roth C, Kang SM, Batty M, Barthélemy M. Structure of Urban Movements: Polycentric Activity and Entangled Hierarchical Flows. PLoS ONE. 2011;6(1):e15923. pmid:21249210
- View Article
- PubMed/NCBI
- Google Scholar
11. Szell M, Sinatra R, Petri G, Thurner S, Latora V. Understanding mobility in a social petri dish. Sci Rep. 2012;2:457. pmid:22708055
- View Article
- PubMed/NCBI
- Google Scholar
12. Song C, Koren T, Wang P, Barabási AL. Modelling the scaling properties of human mobility. Nat Phys. 2010;6:818–823.
- View Article
- Google Scholar
13. Song C, Qu Z, Blumm N, Barabási AL. Limits of predictability in human mobility. Science. 2010;327(5968):1018–1021. pmid:20167789
- View Article
- PubMed/NCBI
- Google Scholar
14. Hufnagel L, Brockmann D, Geisel T. Forecast and control of epidemics in a globalized worlds. Proc. Natl. Acad. Sci. U.S.A. 2004;101(42):15124–15129. pmid:15477600
- View Article
- PubMed/NCBI
- Google Scholar
15. Belik V, Geisel T, Brockmann D. Natural human mobility patterns and spatial spread of infectious diseases. Phys Rev X. 2011;1:011001.
- View Article
- Google Scholar
16. Colizza V, Barrat A, Barthélemy M, Vespignani A. The role of the airline transportation network in the prediction and predictability of global epidemics. Proc. Natl. Acad. Sci. U.S.A. 2006;103(7):2015–2020. pmid:16461461
- View Article
- PubMed/NCBI
- Google Scholar
17. Balcan D, Colizza V, Goncalves B, Hu H, Ramasco JJ, Vespignani A. Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. U.S.A. 2009;106(51):21484–21489. pmid:20018697
- View Article
- PubMed/NCBI
- Google Scholar
18. Batty M. The New Science of Cities. MIT Press; 2013.
19. Zipf GK. The P 1 P 2/D hypothesis: on the intercity movement of persons. American sociological review. 1946;11(6):677–686.
- View Article
- Google Scholar
20. Simini F, González MC, Maritan A, Barabási AL. A universal model for mobility and migration patterns. Nature. 2012;484(7392):96–100. pmid:22367540
- View Article
- PubMed/NCBI
- Google Scholar
21. Stouffer SA. Intervening Opportunities: A Theory Relating Mobility and Distance. American Sociological Review. 1940;5(6):845–867.
- View Article
- Google Scholar
22. Clauset A, Shalizi CR, Newman MEJ. Power-Law Distributions in Empirical Data. SIAM Rev Soc Ind Appl Math. 2009;51(4):661–703.
- View Article
- Google Scholar
23. Sagarra O. Statistical Complex Analysis of Taxi Mobility in San Francisco. M. Sc. Thesis. Universitat Politècnica de Catalunya; 2011.
24. Santi P, Resta G, Szell M, Sobolevsky S, Strogatz SH, Ratti C. Quantifying the Benefits of Vehicle Pooling with Shareability Networks. Proc. Natl. Acad. Sci. U.S.A. 2014;111(37):13290–13294. pmid:25197046
- View Article
- PubMed/NCBI
- Google Scholar
25. Piorkowski M, Sarafijanovic-Djukic N, Grossglauser M. A Parsimonious Model of Mobile Partitioned Networks. Int Conf Commun Syst Netw (COMSNET) IEEE. 2009.1–10.
26. Spieser K, Treleaven K, Zhang R, Frazzoli E, Morton D, Pavone M. Towards a Systematic Approach to the Design and Evaluation of Automated Mobility-on-Demand Systems: a Case Study in Singapore. Road Vehicle Automation. 2014;229–245.
27. Wilson A. A statistical theory of spatial distribution models. Transportation Research. 1967;1(3):253–269.
- View Article
- Google Scholar
28. Erlander S, Stewart NF. The gravity model in transportation analysis: theory and extensions. VSP Utrecht; 1990.
29. Sagarra O, Pérez-Vicente C, Díaz Guilera A. Statistical mechanics of multi-edge networks. Phys Rev E. 2013;88:062806;1–14.
- View Article
- Google Scholar
30. Sagarra O, Font-Clos F, Pérez-Vicente CJ, Díaz-Guilera A. The configuration multi-edge model: Assessing the effect of fixing node strengths on weighted network magnitudes. Europhys Lett. 2014;107(3):38002.
- View Article
- Google Scholar
31. Peng C, Jin X, Wong K, Shi M, Lio P. Collective Human Mobility Pattern from Taxi Trips in Urban Area. PLoS ONE. 2012;7(4):e34487. pmid:22529917
- View Article
- PubMed/NCBI
- Google Scholar
32. Yang Y, Herrera C, Eagle N, González MC. Limits of Predictability in Commuting Flows in the Absence of Data for Calibration. Sci Rep. 2014;4:5662. pmid:25012599
- View Article
- PubMed/NCBI
- Google Scholar
33. Bianconi G, Pin P, Marsili M. Assessing the relevance of node features for network structure. Proc. Natl. Acad. Sci. U.S.A. 2009;106(28):11433–11438. pmid:19571013
- View Article
- PubMed/NCBI
- Google Scholar
34. Lenormand M, Picornell M, Cantú-Ros OG, Tugores A, Louail T, Herranz R, Barthelemy M, Frías-Martinez E, Ramasco J.J. Cross-checking different sources of mobility information. PLoS ONE. 2014;9(8):e105184. pmid:25133549
- View Article
- PubMed/NCBI
- Google Scholar
35. Lenormand M, Huet S, Gargiulo F, Deffuant G. A universal model of commuting networks. PLoS ONE. 2012;7(10):e45985. pmid:23049691
- View Article
- PubMed/NCBI
- Google Scholar
36. Robillard P. Estimating the OD matrix from observed link volumes. Transportation Research. 1975;9(2):123–128.
- View Article
- Google Scholar
37. Sagarra O. ODME: Origin Destination Multi-Edge network package; 2014. Available from: https://github.com/osagarra/ODME_lite.

[ref1] 1. González MC, Hidalgo CA, Barabási AL. Understanding individual human mobility patterns. Nature. 2008;453(7196):779–782. pmid:18528393
View Article
PubMed/NCBI
Google Scholar

[2] View Article

[3] PubMed/NCBI

[4] Google Scholar

[ref2] 2. Blondel VD, Decuyper A, Krings G. A survey of results on mobile phone datasets analysis. arXiv preprint arXiv:1502.03406. 2015.

[ref3] 3. Bazzani A, Giorgini B, Rambaldi S, Gallotti R, Giovannini L. Statistical laws in urban mobility from microscopic GPS data in the area of Florence. J Stat Mech. 2010;2010:P05001
View Article
Google Scholar

[7] View Article

[8] Google Scholar

[ref4] 4. Brockmann D, Hufnagel L, Geisel T. The scaling laws of human travel. Nature. 2006;439(7075):462–465. pmid:16437114
View Article
PubMed/NCBI
Google Scholar

[10] View Article

[11] PubMed/NCBI

[12] Google Scholar

[ref5] 5. Thiemann C, Theis F, Grady D, Brune R, Dirk Brockmann D. The Structure of Borders in a Small World. PLoS ONE. 2010;5:e15422. pmid:21124970
View Article
PubMed/NCBI
Google Scholar

[14] View Article

[15] PubMed/NCBI

[16] Google Scholar

[ref6] 6. Scellato S, Noulas A, Lambiotte R, Mascolo C. Socio-spatial Properties of Online Location-based Social Networks. Proc Int AAAI Conf Weblogs Soc Media. 2011;11.
View Article
Google Scholar

[18] View Article

[19] Google Scholar

[ref7] 7. Scellato S, Musolesi M, Mascolo C, Latora V, Campbell A. NextPlace: A Spatio-Temporal Prediction Framework for Pervasive Systems. Pervasive’11 Proceedings of the 9th international conference on Pervasive computing. 2011;152–169.

[ref8] 8. Barthélemy M. Spatial Networks. Phys Rep. 2010;499:1–101.
View Article
Google Scholar

[22] View Article

[23] Google Scholar

[ref9] 9. Cattuto C, Van den Broeck W, Barrat A, Colizza V, Pinton JF, Vespignani A. Dynamics of person-to-person interactions from distributed RFID sensor networks. PloS ONE. 2010;5(7):e11596. pmid:20657651
View Article
PubMed/NCBI
Google Scholar

[25] View Article

[26] PubMed/NCBI

[27] Google Scholar

[ref10] 10. Roth C, Kang SM, Batty M, Barthélemy M. Structure of Urban Movements: Polycentric Activity and Entangled Hierarchical Flows. PLoS ONE. 2011;6(1):e15923. pmid:21249210
View Article
PubMed/NCBI
Google Scholar

[29] View Article

[30] PubMed/NCBI

[31] Google Scholar

[ref11] 11. Szell M, Sinatra R, Petri G, Thurner S, Latora V. Understanding mobility in a social petri dish. Sci Rep. 2012;2:457. pmid:22708055
View Article
PubMed/NCBI
Google Scholar

[33] View Article

[34] PubMed/NCBI

[35] Google Scholar

[ref12] 12. Song C, Koren T, Wang P, Barabási AL. Modelling the scaling properties of human mobility. Nat Phys. 2010;6:818–823.
View Article
Google Scholar

[37] View Article

[38] Google Scholar

[ref13] 13. Song C, Qu Z, Blumm N, Barabási AL. Limits of predictability in human mobility. Science. 2010;327(5968):1018–1021. pmid:20167789
View Article
PubMed/NCBI
Google Scholar

[40] View Article

[41] PubMed/NCBI

[42] Google Scholar

[ref14] 14. Hufnagel L, Brockmann D, Geisel T. Forecast and control of epidemics in a globalized worlds. Proc. Natl. Acad. Sci. U.S.A. 2004;101(42):15124–15129. pmid:15477600
View Article
PubMed/NCBI
Google Scholar

[44] View Article

[45] PubMed/NCBI

[46] Google Scholar

[ref15] 15. Belik V, Geisel T, Brockmann D. Natural human mobility patterns and spatial spread of infectious diseases. Phys Rev X. 2011;1:011001.
View Article
Google Scholar

[48] View Article

[49] Google Scholar

[ref16] 16. Colizza V, Barrat A, Barthélemy M, Vespignani A. The role of the airline transportation network in the prediction and predictability of global epidemics. Proc. Natl. Acad. Sci. U.S.A. 2006;103(7):2015–2020. pmid:16461461
View Article
PubMed/NCBI
Google Scholar

[51] View Article

[52] PubMed/NCBI

[53] Google Scholar

[ref17] 17. Balcan D, Colizza V, Goncalves B, Hu H, Ramasco JJ, Vespignani A. Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. U.S.A. 2009;106(51):21484–21489. pmid:20018697
View Article
PubMed/NCBI
Google Scholar

[55] View Article

[56] PubMed/NCBI

[57] Google Scholar

[ref18] 18. Batty M. The New Science of Cities. MIT Press; 2013.

[ref19] 19. Zipf GK. The P 1 P 2/D hypothesis: on the intercity movement of persons. American sociological review. 1946;11(6):677–686.
View Article
Google Scholar

[60] View Article

[61] Google Scholar

[ref20] 20. Simini F, González MC, Maritan A, Barabási AL. A universal model for mobility and migration patterns. Nature. 2012;484(7392):96–100. pmid:22367540
View Article
PubMed/NCBI
Google Scholar

[63] View Article

[64] PubMed/NCBI

[65] Google Scholar

[ref21] 21. Stouffer SA. Intervening Opportunities: A Theory Relating Mobility and Distance. American Sociological Review. 1940;5(6):845–867.
View Article
Google Scholar

[67] View Article

[68] Google Scholar

[ref22] 22. Clauset A, Shalizi CR, Newman MEJ. Power-Law Distributions in Empirical Data. SIAM Rev Soc Ind Appl Math. 2009;51(4):661–703.
View Article
Google Scholar

[70] View Article

[71] Google Scholar

[ref23] 23. Sagarra O. Statistical Complex Analysis of Taxi Mobility in San Francisco. M. Sc. Thesis. Universitat Politècnica de Catalunya; 2011.

[ref24] 24. Santi P, Resta G, Szell M, Sobolevsky S, Strogatz SH, Ratti C. Quantifying the Benefits of Vehicle Pooling with Shareability Networks. Proc. Natl. Acad. Sci. U.S.A. 2014;111(37):13290–13294. pmid:25197046
View Article
PubMed/NCBI
Google Scholar

[74] View Article

[75] PubMed/NCBI

[76] Google Scholar

[ref25] 25. Piorkowski M, Sarafijanovic-Djukic N, Grossglauser M. A Parsimonious Model of Mobile Partitioned Networks. Int Conf Commun Syst Netw (COMSNET) IEEE. 2009.1–10.

[ref26] 26. Spieser K, Treleaven K, Zhang R, Frazzoli E, Morton D, Pavone M. Towards a Systematic Approach to the Design and Evaluation of Automated Mobility-on-Demand Systems: a Case Study in Singapore. Road Vehicle Automation. 2014;229–245.

[ref27] 27. Wilson A. A statistical theory of spatial distribution models. Transportation Research. 1967;1(3):253–269.
View Article
Google Scholar

[80] View Article

[81] Google Scholar

[ref28] 28. Erlander S, Stewart NF. The gravity model in transportation analysis: theory and extensions. VSP Utrecht; 1990.

[ref29] 29. Sagarra O, Pérez-Vicente C, Díaz Guilera A. Statistical mechanics of multi-edge networks. Phys Rev E. 2013;88:062806;1–14.
View Article
Google Scholar

[84] View Article

[85] Google Scholar

[ref30] 30. Sagarra O, Font-Clos F, Pérez-Vicente CJ, Díaz-Guilera A. The configuration multi-edge model: Assessing the effect of fixing node strengths on weighted network magnitudes. Europhys Lett. 2014;107(3):38002.
View Article
Google Scholar

[87] View Article

[88] Google Scholar

[ref31] 31. Peng C, Jin X, Wong K, Shi M, Lio P. Collective Human Mobility Pattern from Taxi Trips in Urban Area. PLoS ONE. 2012;7(4):e34487. pmid:22529917
View Article
PubMed/NCBI
Google Scholar

[90] View Article

[91] PubMed/NCBI

[92] Google Scholar

[ref32] 32. Yang Y, Herrera C, Eagle N, González MC. Limits of Predictability in Commuting Flows in the Absence of Data for Calibration. Sci Rep. 2014;4:5662. pmid:25012599
View Article
PubMed/NCBI
Google Scholar

[94] View Article

[95] PubMed/NCBI

[96] Google Scholar

[ref33] 33. Bianconi G, Pin P, Marsili M. Assessing the relevance of node features for network structure. Proc. Natl. Acad. Sci. U.S.A. 2009;106(28):11433–11438. pmid:19571013
View Article
PubMed/NCBI
Google Scholar

[98] View Article

[99] PubMed/NCBI

[100] Google Scholar

[ref34] 34. Lenormand M, Picornell M, Cantú-Ros OG, Tugores A, Louail T, Herranz R, Barthelemy M, Frías-Martinez E, Ramasco J.J. Cross-checking different sources of mobility information. PLoS ONE. 2014;9(8):e105184. pmid:25133549
View Article
PubMed/NCBI
Google Scholar

[102] View Article

[103] PubMed/NCBI

[104] Google Scholar

[ref35] 35. Lenormand M, Huet S, Gargiulo F, Deffuant G. A universal model of commuting networks. PLoS ONE. 2012;7(10):e45985. pmid:23049691
View Article
PubMed/NCBI
Google Scholar

[106] View Article

[107] PubMed/NCBI

[108] Google Scholar

[ref36] 36. Robillard P. Estimating the OD matrix from observed link volumes. Transportation Research. 1975;9(2):123–128.
View Article
Google Scholar

[110] View Article

[111] Google Scholar

[ref37] 37. Sagarra O. ODME: Origin Destination Multi-Edge network package; 2014. Available from: https://github.com/osagarra/ODME_lite.

Figures

Abstract

Introduction

Results

Stability of temporal urban-mobility patterns

The anatomy of urban flows

A flexible model to reproduce human mobility

Supersampling urban trips

Assessing the quality of the supersampling methodology

Discussion

Materials and Methods

Data set

Null model: The Multi-Edge Configuration model

Indicators for the quality of the reconstruction

Distance based measures.

Network measures.

Information values.

Simulations

Acknowledgments

Author Contributions

References