Validate against census flows

Question. Do the home-to-work flows that I measure from trajectories agree with an official table of flows, and which distance decay do they imply?

You will get home and work places from trajectories, expansion weights for a biased sample, the agreement of weighted and unweighted flows with a reference table (common part of commuters, R-squared, Poisson deviance), and a gravity model with a confidence interval for its distance exponent. The last section repeats the gravity and radiation comparison on real check-ins.

Prerequisites. Measure individual mobility and Flows between places.

Data. First a synthetic region of 36 zones with a known answer: a reference ("census") table of commuting flows generated by a gravity model with distance exponent -1.8, and a sample of 27,000 people from it, with a participation rate that grows with the zone number (a biased sample). Official tables such as INSEE's mobility files or the US commuting flows have the same shape (an origin, a destination and a count) and can replace the synthetic table. Then Brightkite check-ins, the first 300,000 rows downloaded with fastmob.data.load_brightkite.

%pip install -q 'fastmob>=0.3.4' pandas statsmodels geopandas shapely

A reference table with a known answer

We lay out 36 zones on a 6 by 6 grid (about 13 km apart), give each a number of resident workers and of jobs, and generate the reference flows with a singly constrained gravity model: the flow from an origin is split among destinations in proportion to their jobs and to distance ** -1.8. Every person in the sample below is drawn from this table, so we can check what the analysis recovers.

import fastmob
import numpy as np
import pandas as pd
from fastmob.core import FlowDataFrame, Locations
from fastmob.measures.evaluation import compare_flow_models
from fastmob.models import Gravity

rng = np.random.default_rng(0)
side = 6
zone = np.arange(side * side)
zones = pd.DataFrame(
    {
        "tile_id": zone,
        "lat": 45.0 + 0.12 * (zone // side),
        "lng": 7.0 + 0.16 * (zone % side),
        "relevance": rng.integers(2_000, 40_000, side * side),  # jobs
        "tot_outflow": rng.integers(20_000, 80_000, side * side),  # resident workers
    }
)
locations = Locations.from_tessellation(
    zones, tile_id_col="tile_id", lat_col="lat", lng_col="lng"
)

truth = Gravity(gravity_type="singly constrained", deterrence_func_args=[-1.8])
census = truth.generate(locations, out_format="flows_sample", random_state=1)
reference = census.to_pandas()
print(
    len(reference),
    "origin-destination pairs,",
    int(reference["flow"].sum()),
    "commuters",
)
reference.nlargest(3, "flow")

1260 origin-destination pairs, 1836919 commuters

The reference table has 1,260 pairs of zones (36 zones, every pair but the diagonal) and about 1.84 million commuters; the largest flows are around 15,000 commuters. The table is on the scale of a whole region, so a sample of a few percent has to be scaled up before it can be compared in absolute terms.

A biased sample of people

A real data set covers a few percent of the population, and not evenly. We draw each commuter of the reference table into the sample with a probability that grows from 0.4% in the first zone to 2.4% in the last, then give each sampled person fixed home and work spots in their zones and six days of records: a night record at home (23:00) and a day record at work (10:00).

participation = 0.004 + 0.02 * zone / zone.max()
drawn = rng.binomial(reference["flow"].astype(int), participation[reference["origin"]])
people = reference.loc[reference.index.repeat(drawn), ["origin", "destination"]]
people = people.reset_index(drop=True).assign(uid=lambda d: d.index)


def spot(zone_of_person):
    return (
        45.0 + 0.12 * (zone_of_person // side) + rng.uniform(-0.03, 0.03, len(people)),
        7.0 + 0.16 * (zone_of_person % side) + rng.uniform(-0.04, 0.04, len(people)),
    )


home_lat, home_lng = spot(people["origin"].to_numpy())
work_lat, work_lng = spot(people["destination"].to_numpy())
monday = pd.Timestamp("2024-03-04")
records = []
for day in range(6):
    for hour, lat, lng in ((23, home_lat, home_lng), (10, work_lat, work_lng)):
        when = monday + pd.Timedelta(days=day, hours=hour)
        records.append(
            pd.DataFrame(
                {"uid": people["uid"], "datetime": when, "lat": lat, "lng": lng}
            )
        )
trajectories = pd.concat(records, ignore_index=True)
print(len(people), "people,", len(trajectories), "records")

27420 people, 329040 records

Home and work from trajectories

home_location is the most visited place at night and work_location the most visited place during working hours on weekdays. We map each place to its zone (the grid makes this a rounding) and count the people for each pair of zones. People whose home and work fall in the same zone are not commuters between zones and are left out, as in the reference table.

home = fastmob.home_location(trajectories)
work = fastmob.work_location(trajectories)


def zone_of(places):
    row = np.rint((places["lat"] - 45.0) / 0.12).astype(int).clip(0, side - 1)
    col = np.rint((places["lng"] - 7.0) / 0.16).astype(int).clip(0, side - 1)
    return (row * side + col).to_numpy()


pairs = home.assign(home=zone_of(home))[["uid", "home"]].merge(
    work.assign(work=zone_of(work))[["uid", "work"]], on="uid"
)
pairs = pairs[pairs["home"] != pairs["work"]]
print(len(pairs), "people with different home and work zones")
pairs.groupby("home").size().describe().round(0)

27420 people with different home and work zones





count      36.0
mean      762.0
std       434.0
min        94.0
25%       364.0
50%       710.0
75%      1167.0
max      1649.0
dtype: float64

The sample has 27,420 people, 762 per origin zone on average, from 94 to 1,649 depending on the zone. The sampling rate grows with the zone number, as designed.

Weighting the sample

A person from a zone with many sampled people stands for fewer commuters than a person from a zone with few. The expansion weight of everyone living in zone h is residents(h) / sampled(h): the number of resident workers in the reference table over the number of sampled people. Summing the weights of the people of each pair of zones gives the estimated flow. As a baseline we use the plain counts, multiplied by one constant so that the total matches (otherwise the flows are on the scale of the sample and cannot be compared in absolute terms).

residents = zones.set_index("tile_id")["tot_outflow"]
sampled = pairs.groupby("home").size()
pairs = pairs.assign(weight=pairs["home"].map(residents / sampled))
constant = reference["flow"].sum() / len(pairs)


def flow_table(weights):
    summed = (
        pairs.assign(flow=weights)
        .groupby(["home", "work"], as_index=False)["flow"]
        .sum()
    )
    return summed.rename(columns={"home": "origin", "work": "destination"})


def as_flows(table):
    return FlowDataFrame(table, locations=locations, tile_id="tile_id")


estimates = {
    "unweighted": as_flows(flow_table(constant)),
    "weighted": as_flows(flow_table(pairs["weight"])),
}
pairs["weight"].describe().round(0)

count    27420.0
mean        67.0
std         31.0
min         43.0
25%         46.0
50%         58.0
75%         70.0
max        264.0
Name: weight, dtype: float64

The weights run from 43 to 264 commuters per person, with a median of 58: a person in a lightly sampled zone stands for many more commuters than one in a heavily sampled zone. The constant used for the unweighted baseline is 67 commuters per person (the average weight).

Agreement with the reference

compare_flow_models scores each estimate against the reference over a common set of origin-destination pairs (every pair of the reference table's locations by default). The common part of commuters (CPC) is the share of commuters that both tables send along the same pair of zones (1 is identical); R-squared and the rank correlation compare the pair flows; the Poisson deviance and the pseudo R-squared treat the flows as counts.

reference_flows = as_flows(reference[["origin", "destination", "flow"]])
scores = compare_flow_models(reference_flows, estimates).to_pandas().set_index("model")
scores[["n_pairs", "cpc", "cpl", "r_squared", "spearman", "pseudo_r_squared"]].round(3)

Weighting raises every score. The CPC goes from 0.83 to 0.92 (92% of the commuters are placed on the same pair of zones as in the reference), the R-squared of the pair flows from 0.86 to 0.98 and the rank correlation from 0.90 to 0.95. The CPL is 0.97 for both, because the sets of pairs with a flow are the same: weighting changes the sizes of the flows, not which pairs have flow. The pseudo R-squared is lower than the R-squared because they answer different questions: the R-squared follows the large flows, while the Poisson deviance counts the error on every pair, including the many small ones.

The distance decay, with a confidence interval

Gravity.fit fits the model by Poisson regression and fit_summary returns each parameter with a robust standard error, a 95% confidence interval and a p-value. We fit the reference table itself (to check the method), the unweighted estimate and the weighted one. The singly constrained model has one free level per origin, so a bias that depends only on the origin zone, like ours, does not bias the distance exponent when the model holds.

rows = []
for name, flows in {"reference": reference_flows, **estimates}.items():
    model = Gravity(gravity_type="singly constrained")
    model.fit(flows)
    fit = model.fit_summary().to_pandas().set_index("parameter")
    stats = model.fit_statistics()
    rows.append(
        {
            "table": name,
            **fit.loc["distance_coefficient", ["estimate", "ci_low", "ci_high"]],
            "dispersion": stats["dispersion"],
            "pseudo_r_squared": stats["pseudo_r_squared"],
        }
    )
pd.DataFrame(rows).set_index("table").round(3)

The reference table gives -1.800 with an interval of -1.803 to -1.798, so the method recovers the exponent that generated it (the dispersion of 0.96 is the Poisson value). The unweighted sample gives -1.797 (interval -1.821 to -1.774) and the weighted one -1.791 (-1.816 to -1.765): both within 0.01 of the truth, with the truth inside both intervals. The weights move the estimate a little (by 0.006) and its precision, which in theory they should not bias; their job is the flow levels, as the scores above show. The dispersions of 69 and 91 are mostly the size of the weights (a person stands for 67 commuters), not behaviour, so do not read them as overdispersion here. The robust intervals are the ones to use.

With a real census table

Our zones are squares of a grid, so a rounding gave the zone of a place. With real zones (communes, counties, census tracts) you need a point-in-polygon join and the same zone identifiers as in the census table. With GeoPandas, sjoin does it; we check it against the rounding on our squares. Real polygons often come in a projected system (French communes in Lambert-93, EPSG:2154), so convert the points or the polygons to the same coordinate system first (to_crs). The census table must have an origin, a destination and a count: for the French census (INSEE, the mobility file MOBPRO) the commune of residence, the commune of work and the weight of each record; check the layout of the file you download. The weights need the number of resident workers per zone, and since the sample leaves out people who work in their home zone, take it from the census flows between different zones (the sum over destinations other than the origin).

import geopandas as gpd
from shapely.geometry import box

squares = gpd.GeoDataFrame(
    {"zone": zone},
    geometry=[
        box(
            7.0 + 0.16 * c - 0.08,
            45.0 + 0.12 * r - 0.06,
            7.0 + 0.16 * c + 0.08,
            45.0 + 0.12 * r + 0.06,
        )
        for r, c in zip(zone // side, zone % side)
    ],
    crs="EPSG:4326",
)
points = gpd.GeoDataFrame(
    home, geometry=gpd.points_from_xy(home["lng"], home["lat"]), crs="EPSG:4326"
)
joined = gpd.sjoin(points, squares, predicate="within")
same = joined["zone"].reindex(points.index).to_numpy() == zone_of(home)
print(f"point-in-polygon and rounding agree for {same.mean():.0%} of the home places")

point-in-polygon and rounding agree for 100% of the home places

Exactly the same pairs follow, so the rest of the analysis is unchanged: the pairs, the weights residents / sampled, compare_flow_models and Gravity.fit.

The same on real trajectories

For real data there is no census table, but the questions are the same: which distance decay do the flows show, and which model reproduces them better? We take the Brightkite check-ins, group them into H3 cells of about 22 km (resolution 4), and keep the 60 busiest cells and the trips between them. The relevance of a cell is its number of check-ins, which is itself made of the trips we model: an independent size such as population would be better, and the destination exponent here is partly circular. trajectory_to_od lists only the pairs that have trips; Gravity.fit adds the other pairs as zeros, which matters (below).

from fastmob.data import load_brightkite

raw = load_brightkite(max_rows=300_000).to_pandas()
raw = raw[~(((raw["lat"] == 0) & (raw["lng"] == 0)) | (raw["lat"].abs() >= 90))]
counts = raw.groupby("uid").size()
checkins = (
    raw[raw["uid"].isin(counts[counts >= 100].index)]
    .sort_values(["uid", "datetime"])
    .reset_index(drop=True)
)


def busiest_cells(frame):
    """The 60 cells with most check-ins, with their number of check-ins as relevance."""
    counted = fastmob.latlng_to_h3(frame, resolution=4).groupby("h3_cell").size()
    return counted.rename("relevance").nlargest(60).rename_axis("tile_id")


def cell_flows(frame, busiest):
    """Observed flows between the busiest cells, for the given check-ins, with their catalogue."""
    trips = fastmob.trajectory_to_od(frame, resolution=4).rename(
        columns={"count": "flow"}
    )
    trips = trips[
        trips["origin"].isin(busiest.index)
        & trips["destination"].isin(busiest.index)
        & (trips["origin"] != trips["destination"])
    ].reset_index(drop=True)
    centres = fastmob.h3_to_latlng(pd.DataFrame({"h3_cell": busiest.index})).rename(
        columns={"h3_cell": "tile_id", "center_lat": "lat", "center_lng": "lng"}
    )
    table = (
        centres.join(busiest, on="tile_id")
        .join(trips.groupby("origin")["flow"].sum().rename("tot_outflow"), on="tile_id")
        .fillna({"tot_outflow": 0})
    )
    catalogue = Locations.from_tessellation(
        table, tile_id_col="tile_id", lat_col="lat", lng_col="lng"
    )
    return FlowDataFrame(trips, locations=catalogue, tile_id="tile_id"), catalogue


observed, real_locations = cell_flows(checkins, busiest_cells(checkins))
print(len(observed.to_pandas()), "pairs with trips out of", 60 * 59, "pairs of cells")

646 pairs with trips out of 3540 pairs of cells

rows = {}
for label, include in (("pairs with trips only", False), ("zero pairs added", True)):
    model = Gravity(gravity_type="singly constrained")
    model.fit(observed, include_zeros=include)
    fit = model.fit_summary().to_pandas().set_index("parameter")
    rows[label] = {
        "flows used": model.fit_statistics()["n_observations"],
        **fit.loc["distance_coefficient", ["estimate", "ci_low", "ci_high"]],
        "dispersion": model.fit_statistics()["dispersion"],
    }
pd.DataFrame(rows).T.round(3)

trajectory_to_od lists only the 646 pairs of cells that have trips. Four of the 60 cells have no outgoing trip, so 56 origins times 59 destinations make 3,304 pairs, and 2,658 of them have no trip. Fitted on the 646 alone, the distance coefficient is -1.21 (interval -1.27 to -1.15). With the other 2,658 pairs added as zeros, which is what Gravity.fit does by default, it is -1.40 (-1.47 to -1.33), and the intervals do not overlap: leaving out the pairs that nobody travelled flattens the decay. The dispersion falls from 15.8 to 6.6.

To compare gravity and radiation fairly we cannot score a model on the flows it was fitted to: gravity has parameters fitted to them and radiation has none. We split the users into two halves, fit the gravity model on the flows of the first half, and score both models against the flows of the second half. Both are generated with the outflows of the second half, so they differ only in how they distribute the trips.

from fastmob.models import Radiation

users = np.sort(checkins["uid"].unique())
first = checkins["uid"].isin(users[::2])
train, test = checkins[first], checkins[~first]
cells_of_first_half = busiest_cells(
    train
)  # chosen, like the relevance, without the second half
train_flows, _ = cell_flows(train, cells_of_first_half)
test_flows, test_locations = cell_flows(test, cells_of_first_half)

gravity = Gravity(gravity_type="singly constrained")
gravity.fit(train_flows)
second = Gravity(gravity_type="singly constrained")
second.fit(test_flows)
pd.concat(
    {
        name: fitted.fit_summary()
        .to_pandas()
        .set_index("parameter")[["estimate", "ci_low", "ci_high"]]
        for name, fitted in (("first half", gravity), ("second half", second))
    }
).round(3)

The two halves of the users give distance coefficients of -1.54 (interval -1.63 to -1.45) and -1.42 (-1.50 to -1.35), and destination exponents of 0.53 and 0.67: similar, but not the same, which is the size of the sampling difference between two groups of 150 people. Both differ from the -1.40 of all the users above, because the cells and their relevance are chosen again for each group.

models = {
    "gravity (fitted on the other half)": gravity.generate(
        test_locations, out_format="flows"
    ),
    "radiation": Radiation().generate(test_locations, out_format="flows"),
}
compare_flow_models(test_flows, models).to_pandas().set_index("model").round(3)

On the second half, the gravity model fitted on the first half reproduces the flows better than the radiation model, which has no parameters to fit: CPC 0.79 against 0.55, R-squared 0.97 against 0.77, pseudo R-squared 0.93 against 0.58 and a Poisson deviance six times smaller; the rank correlations are close (0.43 and 0.40). Both models are generated with the outflows of the second half, so the comparison is about how the trips are distributed, and the cells and relevance come from the first half only. The radiation model has a higher CPL (0.45 against 0.24), but a model that gives every pair some flow has a CPL that is fixed by the number of observed pairs, so do not read it. All scores are on the same 3,540 pairs (every pair of the 60 cells). These are 60 cells of 22 km in a sample of check-ins, so the numbers describe this data set, not human mobility in general.

Interpretation and limits

• The reference table here is synthetic, which is what lets us check that the method recovers the truth. With an official table the check is only as good as that table: commuting tables count where people work or live by registration, and sampled trajectories can include second homes, remote work or night-shift work.
• The expansion weight assumes that, within a home zone, the people sampled are like the other residents. A bias that depends on who works where (not only on where people live) is not corrected, and shows as a lower CPC.
• Home and work are inferred. Change start_night, end_night, start_work and end_work, or the number of nights required, and see whether the flows change.
• CPC and the other scores depend on the size of the zones: finer zones give lower values. Report the zone size with the score.
• The Poisson regression assumes counts. Real flows are overdispersed (the dispersion above is far above 1), so use the robust intervals that fit_summary reports, and treat p-values as approximate.
• A gravity fit on a few dozen cells describes those cells. The distance exponent depends on the range of distances and on the zone size.

Next steps

• Fit and compare the distributions of jump lengths with intervals and alternatives: Compare to mobility laws.
• How the fits are estimated and what to report: How the fits work.
• Score a flow model with more metrics: Fit and compare distributions.

Related API

home_location · work_location · Gravity.fit · Gravity.fit_summary · Gravity.fit_statistics · Radiation.generate · compare_flow_models