How to estimate Poisson origin–destination flow models¶

This guide fits SARPoissonFlowSeparable to origin–destination counts and decides between Poisson and negative binomial observation noise for the same spatial structure.

Use it when your flows are counts, you want a correctly specified count likelihood rather than a log-transformed Gaussian, and you need to know whether the extra dispersion parameter of a negative binomial is buying you anything.

The Poisson flow models sample by auxiliary-mixture Gibbs (Frühwirth-Schnatter & Wagner 2006): Poisson inter-arrival times are augmented and \(-\log\tau\) is approximated by a ten-component normal mixture, which leaves a conditionally Gaussian model that the existing SAR machinery — the same \(\rho\) step, the same sparse solves, the same log-determinant — already serves.

Simulate flow counts¶

generate_poisson_flow_data_separable draws counts from \(y_{ij}\sim\text{Poisson}(\mu_{ij})\) with \(\log\boldsymbol\mu = A(\rho_d,\rho_o,\rho_w)^{-1}X\beta\) and the separable restriction \(\rho_w = -\rho_d\rho_o\). With \(n = 25\) regions there are \(N = n^2 = 625\) flows.

import time

import arviz as az
import numpy as np
import pandas as pd

from neighbayes.dgp.flows import generate_poisson_flow_data_separable
from neighbayes.models.flow import SARNegBinFlowSeparable, SARPoissonFlowSeparable

RHO_D, RHO_O = 0.35, 0.25

data = generate_poisson_flow_data_separable(n=25, rho_d=RHO_D, rho_o=RHO_O, seed=7)
y, X, G = data["y_vec"], data["X"], data["G"]

print(f"flows      : {y.size}")
print(f"mean count : {y.mean():.2f}")
print(f"variance   : {y.var():.2f}")
print(f"zeros      : {(y == 0).mean():.1%}")

flows      : 625
mean count : 8.12
variance   : 51.81
zeros      : 4.5%

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Note the variance-to-mean ratio. These counts are conditionally Poisson by construction, yet marginally they look badly overdispersed, because the spatial filter \(A^{-1}\) spreads \(\mu_{ij}\) over orders of magnitude across cells. Marginal overdispersion is not evidence that you need a negative binomial — in flow data it is the ordinary signature of a gravity mean with spatial feedback. The question has to be settled by fitting both.

Fit the Poisson flow model¶

The interface is the shared flow-model interface: pass the vectorized flow counts, the \(N\times p\) origin–destination design matrix, and the regional weights graph on \(n\) units. Poisson flow models are Gibbs-only — fit(sampler="nuts") raises.

model_pois = SARPoissonFlowSeparable(y, X, G, col_names=data["col_names"])

start = time.perf_counter()
idata_pois = model_pois.fit(
    draws=1000, tune=500, chains=2, random_seed=11, progressbar=False
)
secs_pois = time.perf_counter() - start

az.summary(idata_pois, var_names=["rho_d", "rho_o", "rho_w"])

rho_w is deterministic here — the separable model pins it to \(-\rho_d\rho_o\) rather than sampling it — so it carries no independent information. It is reported because the spatial-effects decomposition needs all three.

Check recovery against the truth:

post = idata_pois.posterior
print(f"rho_d: {float(post['rho_d'].mean()):.3f}  (true {RHO_D})")
print(f"rho_o: {float(post['rho_o'].mean()):.3f}  (true {RHO_O})")

rho_d: 0.311  (true 0.35)
rho_o: 0.226  (true 0.25)

Read the spillovers¶

spatial_effects() gives the LeSage–Thomas-Agnan decomposition of each covariate into origin, destination, intra-regional and network components.

model_pois.spatial_effects().round(4)

Decide between Poisson and negative binomial¶

Fit SARNegBinFlowSeparable to the same counts, the same design, the same weights. Only the observation model differs — the negative binomial adds a free dispersion parameter alpha.

model_nb = SARNegBinFlowSeparable(y, X, G, col_names=data["col_names"])

start = time.perf_counter()
idata_nb = model_nb.fit(
    draws=1000, tune=500, chains=2, random_seed=11, progressbar=False
)
secs_nb = time.perf_counter() - start

az.summary(idata_nb, var_names=["rho_d", "rho_o", "alpha"])

Look at alpha before anything else. Under the negative binomial, \(\operatorname{Var}(y) = \mu + \mu^2/\alpha\), so \(\alpha \to \infty\) is the Poisson. A posterior mean in the hundreds with a standard deviation of the same order is the model telling you it cannot locate the dispersion parameter, because there is no conditional overdispersion to locate. That is the first piece of evidence for the simpler likelihood.

Its converse is worth knowing too: as \(\alpha\) grows, the Pólya–Gamma augmentation the negative binomial sampler relies on degenerates — the working precision \(E[\omega]\) diverges while the marginal Fisher information stays at \(\mu\), and effective sample size collapses. Fixing \(\alpha\) at a large value to “approximate” a Poisson is therefore not a shortcut; it is the one regime where that sampler fails. The auxiliary-mixture scheme exists to reach the Poisson directly instead.

Do the two posteriors agree?¶

Two independently implemented samplers landing on the same \(\rho\) posterior is mutual validation. Disagreement beyond Monte Carlo error is a specification signal, not a rounding difference.

rows = []
for name in ("rho_d", "rho_o"):
    rows.append(
        {
            "parameter": name,
            "truth": {"rho_d": RHO_D, "rho_o": RHO_O}[name],
            "poisson": float(idata_pois.posterior[name].mean()),
            "negbin": float(idata_nb.posterior[name].mean()),
            "poisson sd": float(idata_pois.posterior[name].std()),
            "negbin sd": float(idata_nb.posterior[name].std()),
        }
    )
pd.DataFrame(rows).set_index("parameter").round(4)

What does each cost?¶

Sampling efficiency is effective sample size per second, not draws per second. The two likelihoods have different per-sweep costs — the auxiliary-mixture scheme draws two latent inter-arrival times per observation, the Pólya–Gamma scheme draws one mixing weight — so equal draws are not equal work.

rows = []
for label, idata, secs in (
    ("Poisson", idata_pois, secs_pois),
    ("NegBin", idata_nb, secs_nb),
):
    ess = az.ess(idata, var_names=["rho_d", "rho_o"])
    rhat = az.rhat(idata, var_names=["rho_d", "rho_o"])
    rows.append(
        {
            "model": label,
            "seconds": secs,
            "ess rho_d": float(ess["rho_d"]),
            "ess rho_o": float(ess["rho_o"]),
            "ess/sec rho_d": float(ess["rho_d"]) / secs,
            "max rhat": float(max(rhat["rho_d"], rhat["rho_o"])),
        }
    )
pd.DataFrame(rows).set_index("model").round(3)

Neither likelihood dominates on speed, and which one comes out ahead shifts with \(n\), with the count level, and with how much of the mass sits near zero — so measure it on your own data rather than carrying a rule of thumb. Efficiency is the tiebreaker, not the criterion. Both samplers here recover the same \(\rho\) posterior with \(\hat R\) at 1.00, so nothing about the science turns on the choice; what should decide it is whether the observation model is right, which is the predictive check below.

Posterior predictive check in count space¶

The honest specification test is whether replicated counts reproduce features of the observed data that the likelihood did not fit directly — here the variance-to-mean ratio and the share of zeros.

y_rep = model_pois.posterior_predictive(random_seed=0)

obs = {"var/mean": y.var() / y.mean(), "zeros": (y == 0).mean(), "max": y.max()}
rep = {
    "var/mean": np.mean(y_rep.var(axis=1) / y_rep.mean(axis=1)),
    "zeros": np.mean((y_rep == 0).mean(axis=1)),
    "max": np.mean(y_rep.max(axis=1)),
}
pd.DataFrame({"observed": obs, "replicated (mean)": rep}).round(3)

If the replicated variance-to-mean ratio and zero share bracket the observed values, the Poisson is doing its job and the negative binomial’s dispersion parameter is fitting noise. If replicated counts are systematically tighter than observed — too little spread, too few zeros — that is genuine excess dispersion and the negative binomial earns its extra parameter.

For real trade and migration flows, Santos Silva & Tenreyro (2011) argue that zero inflation is usually the wrong fix for excess zeros: the zeros in gravity data are ordinary Poisson zeros from small \(\mu_{ij}\), not a separate structural process. Running the comparison above on your own data is the way to check that claim rather than assume it.

When you want all three \(\rho\)¶

SARPoissonFlow exposes the unrestricted three-parameter model, with \(\rho_w\) free rather than pinned to \(-\rho_d\rho_o\).

See also¶

• How to estimate origin–destination flow models — the Gaussian and negative binomial flow families, and the effects decomposition in full

• How to estimate panel flow models — the same models over repeated periods

• Supported Models — the full model catalogue