Poisson regression, negative binomial and PPML gravity models in Python

Counts such as doctor visits, patents or accidents are nonnegative whole numbers, often with many zeros and a long right tail. Trade flows between countries are amounts rather than counts, but they share the problem: they cannot be negative and many are exactly zero. A linear regression can predict negative values. A regression of log(y) has to drop the zeros, and under heteroskedasticity it estimates the wrong parameters (Santos Silva and Tenreyro 2006). Poisson regression and its relatives model the logarithm of the conditional mean instead.

What Poisson, negative binomial and PPML estimate

All three models share the same exponential mean:

E[yᵢ | xᵢ] = μᵢ = exp(xᵢ′β),   Var(yᵢ | xᵢ) = μᵢ (Poisson),  μᵢ(1 + αμᵢ) (NB2),  μᵢ(1 + δ) (NB1)

A coefficient is a semi-elasticity: one more unit of x changes the expected count by roughly 100·β percent, and exp(β) is the incidence-rate ratio. When the regressor is in logs, such as log distance, its coefficient is an elasticity.

Poisson regression is maximum likelihood under the Poisson distribution, whose variance equals its mean. Real counts are usually overdispersed: the variance is larger. The Poisson estimates nevertheless stay consistent whenever the conditional mean is correctly specified, whatever the distribution of y (Gourieroux, Monfort and Trognon 1984). What overdispersion breaks are the conventional standard errors, which become too small; robust or cluster-robust standard errors repair them. Cameron and Trivedi (1990) give regression-based tests for overdispersion.

The negative binomial model adds a dispersion parameter. NB2, the usual form, has variance μ(1 + αμ); NB1 has μ(1 + δ). When its variance form is right, the negative binomial is more efficient than Poisson. With an estimated dispersion it does not, in general, share the Poisson model's robustness to a misspecified distribution, so Poisson with robust standard errors is the safer choice when only the mean matters.

PPML (Poisson pseudo-maximum likelihood) applies the same estimator to any nonnegative outcome, counts or not. Santos Silva and Tenreyro (2006) showed that it is the natural estimator for gravity equations of trade: it keeps the zero flows and stays consistent under heteroskedasticity, where OLS on log trade does not. Gravity models need many fixed effects, for example exporter and importer effects or exporter-year, importer-year and pair effects. ppmlhdfe absorbs them without forming dummy variables.

In OpenEconometrics

oe.poisson(data=df, y="visits", x=["income", "chronic"]) corresponds to the Stata command poisson visits income chronic:

• exposure="pyears" adds ln(exposure) to the index with coefficient 1, so that the model is for a rate; offset= adds a column directly. categorical=[...] expands columns into treatment dummies.
• covariance="robust" gives sandwich standard errors and cluster= (one or two columns) cluster-robust ones; the default "nonrobust" uses the inverse information.
• tests["gof_deviance"] and tests["gof_pearson"] are the goodness-of-fit tests of Stata's estat gof. tests["model"] is a likelihood-ratio test against the constant-only model under nonrobust, and a Wald test otherwise.

oe.nbreg takes the same arguments plus dispersion="mean" (NB2, the default) or "constant" (NB1, Stata's nbreg ..., dispersion(constant)). The ancillary parameter is reported as the term /lnalpha (or /lndelta). extra["alpha"] gives α with a delta-method standard error and an interval obtained by exponentiating the interval for ln α. tests["alpha"] is the likelihood-ratio test of α = 0 against Poisson, with the p-value halved because the null lies on the boundary; it is omitted under robust and cluster.

oe.ppmlhdfe(data=flows, y="trade", x=["ldist", "rta"], absorb=["exporter", "importer"]) corresponds to the community Stata command ppmlhdfe trade ldist rta, absorb(exporter importer). Any number of fixed-effect dimensions can be absorbed. The default covariance is robust, and cluster= takes one or two columns. Observations in a fixed-effect level whose outcomes are all zero are dropped, because that effect diverges, and so are singletons (drop_singletons=True); both counts are reported. No intercept is shown, since the effects absorb it.

After poisson and nbreg, oe.margins(fit, ["income"], data=df) gives average effects on the expected count. The GLM and count guide documents all three estimators; zero-inflated, truncated and hurdle count models are in the count-model guide.

Supported options

Options of oe.nbreg

Options of oe.ppmlhdfe

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script has two parts. First, it simulates doctor visits for 2,000 people from a negative binomial (NB2) model with mean exp(0.5 − 0.2·income + 0.6·chronic) and α = 0.8, so the data are overdispersed; income is standard normal and half the people have a chronic condition. Second, it simulates trade among 40 countries, all 1,560 ordered pairs, with exporter and importer effects, a distance elasticity of −1.0 and a regional-trade-agreement effect of 0.5. Gamma noise with mean 1 makes the flows heteroskedastic, and many small flows are zero.

import numpy as np
import pandas as pd
import openecon as oe


def doctor_visits(n=2000, seed=5):
    """Overdispersed counts: negative binomial (NB2) with alpha = 0.8."""
    rng = np.random.default_rng(seed)
    income = rng.normal(size=n)
    chronic = rng.integers(0, 2, size=n)
    mu = np.exp(0.5 - 0.2 * income + 0.6 * chronic)
    visits = rng.poisson(mu * rng.gamma(shape=1 / 0.8, scale=0.8, size=n))
    return pd.DataFrame({"visits": visits, "income": income, "chronic": chronic})


def trade_flows(countries=40, seed=6):
    """Gravity: exporter and importer effects, distance elasticity -1.0, RTA effect 0.5."""
    rng = np.random.default_rng(seed)
    xy = rng.uniform(0, 10, size=(countries, 2))
    exporter_fe, importer_fe = rng.normal(size=countries), rng.normal(size=countries)
    i, j = np.where(~np.eye(countries, dtype=bool))
    ldist = np.log(1 + np.hypot(*(xy[i] - xy[j]).T))   # distance in thousand km, at least 1
    rta = (rng.uniform(size=i.size) < 0.2).astype(int)
    mean = np.exp(2.0 + exporter_fe[i] + importer_fe[j] - 1.0 * ldist + 0.5 * rta)
    trade = rng.poisson(mean * rng.gamma(shape=2.0, scale=0.5, size=i.size))
    return pd.DataFrame({"trade": trade, "ldist": ldist, "rta": rta, "exporter": i, "importer": j})


if __name__ == "__main__":
    df = doctor_visits()
    x = ["income", "chronic"]
    poisson = oe.poisson(data=df, y="visits", x=x)
    robust = oe.poisson(data=df, y="visits", x=x, covariance="robust")
    nb = oe.nbreg(data=df, y="visits", x=x)
    truth = {"Intercept": 0.5, "income": -0.2, "chronic": 0.6}
    rows = [(term, truth[term], p.estimate, p.std_error, r.std_error, b.estimate, b.std_error)
            for term, p, r, b in zip(truth, poisson.coefficients, robust.coefficients, nb.coefficients)]
    columns = ["term", "true", "poisson", "se", "robust se", "nbreg", "nbreg se"]
    print(pd.DataFrame(rows, columns=columns).round(4).to_string(index=False))
    gof = poisson.tests["gof_pearson"]
    print(f"Poisson Pearson goodness of fit: chi2({gof['df']}) = {gof['statistic']:.1f}, p = {gof['p_value']:.2g}")
    alpha, lr = nb.extra["alpha"], nb.tests["alpha"]
    print(f"nbreg alpha = {alpha['estimate']:.3f} (95% CI {alpha['ci_low']:.3f} to {alpha['ci_high']:.3f}); "
          f"LR test of alpha = 0: chibar2 = {lr['statistic']:.1f}, p = {lr['p_value']:.2g}")

    flows = trade_flows()
    print(f"\n{len(flows)} country pairs, {(flows.trade == 0).mean():.0%} with zero trade")
    ppml = oe.ppmlhdfe(data=flows, y="trade", x=["ldist", "rta"], absorb=["exporter", "importer"])
    print(ppml.summary())
    positive = flows[flows.trade > 0].assign(ltrade=lambda d: np.log(d.trade))
    ols = oe.reghdfe(data=positive, y="ltrade", x=["ldist", "rta"], absorb=["exporter", "importer"])
    print("OLS of log trade on positive flows:", {c.term: round(c.estimate, 4) for c in ols.coefficients})

With openecon 0.3.20a1 from PyPI the script prints:

     term  true  poisson     se  robust se   nbreg  nbreg se
Intercept   0.5   0.4580 0.0255     0.0394  0.4531    0.0397
   income  -0.2  -0.1876 0.0145     0.0272 -0.1980    0.0267
  chronic   0.6   0.7099 0.0307     0.0528  0.7169    0.0523
Poisson Pearson goodness of fit: chi2(1997) = 6520.4, p = 0
nbreg alpha = 0.880 (95% CI 0.794 to 0.976); LR test of alpha = 0: chibar2 = 1772.3, p = 0

1560 country pairs, 39% with zero trade
Poisson pseudo-likelihood regression with fixed effects — trade
Observations: 1560  |  Covariance: robust  |  Confidence: 95%

Term   Estimate  Std. error         z     P>|stat|  CI lower   CI upper
-----  --------  ----------  --------  -----------  --------  ---------
ldist  -0.95109   0.0971357  -9.79135   1.2265e-22  -1.14147  -0.760707
rta    0.539225    0.106518   5.06231  4.14212e-07  0.330454   0.747996

pseudo_r_squared: 0.656961  |  log_likelihood: -3507.98  |  deviance: 4148.59  |  df_resid: 1479
df_absorbed: 79  |  iterations: 7  |  n_separated_dropped: 0  |  n_singletons_dropped: 0
Wald chi2 test of the slopes: chi2(2) = 118.148, p = 2.211e-26
OLS of log trade on positive flows: {'ldist': -0.655, 'rta': 0.3638}

Reading the output

• The counts are clearly overdispersed. The Pearson statistic of 6520.4 on 1,997 degrees of freedom is more than three times its expected value, the likelihood-ratio test of α = 0 gives 1772.3, and the estimated α of 0.880 has an interval from 0.794 to 0.976 that contains the true 0.8.
• Poisson's conventional standard errors are too small. For the chronic condition the conventional standard error is 0.0307, but the robust one is 0.0528, 1.7 times larger and almost the same as the negative binomial's 0.0523. Here the negative binomial buys little precision over Poisson with robust errors.
• Both models give similar coefficients, and both draw the chronic effect high. The estimates of 0.7099 and 0.7169 lie about two standard errors above the true 0.6. The robust Poisson interval, from 0.606 to 0.813, just contains it; the negative binomial interval, from 0.614 to 0.819, just misses. That is sampling error shared by two estimators of the same mean, not bias.
• PPML recovers the gravity coefficients. The distance elasticity is −0.951 with an interval from −1.141 to −0.761 around the true −1.0. The agreement coefficient of 0.539 means exp(0.539) − 1 ≈ 71% more trade (true 65%). The 79 absorbed degrees of freedom are 40 exporter plus 40 importer effects, less one redundant level.
• OLS on log trade is biased towards zero. Dropping the 39% of zero flows and regressing log trade gives a distance elasticity of −0.655 and an agreement effect of 0.364. Dropping zeros keeps the small flows that happened to draw high noise, and when the variance depends on the regressors the mean of log trade is not the log of mean trade (Santos Silva and Tenreyro 2006).

Limitations

• Cluster-robust covariance for oe.poisson allows at most 2 clustering dimensions.
• Cluster-robust covariance for oe.nbreg allows at most 2 clustering dimensions.
• Cluster-robust covariance for oe.ppmlhdfe allows at most 2 clustering dimensions.
• oe.ppmlhdfe does not search for observations separated by the regressors themselves. Such separation raises separation_detected naming the regressor; drop it or the separated observations yourself.
• oe.ppmlhdfe absorbs columns, not interaction syntax: build a column such as exporter_year for an exporter-by-year effect.
• oe.ppmlhdfe uses reghdfe's finite-sample factors, N/(N − K) for robust and G/(G − 1)·(N − 1)/(N − K) for cluster covariance, which the guide lists as an unconfirmed convention.
• oe.nbreg raises boundary_solution when the data are not overdispersed, rather than reporting α near zero; use oe.poisson, or oe.glm(family="nbinomial", dispersion=...) to fix α.
• Negative binomial likelihoods are evaluated only where y + 1/α + 1 ≤ 10⁷ for NB2 (μ/δ for NB1); a larger individual count raises precision_unsupported. This limits count size, not the number of observations.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Cameron, A. C. and Trivedi, P. K. (1990). Regression-based tests for overdispersion in the Poisson model. Journal of Econometrics, 46(3), 347–364. doi:10.1016/0304-4076(90)90014-K
• Gourieroux, C., Monfort, A. and Trognon, A. (1984). Pseudo maximum likelihood methods: Applications to Poisson models. Econometrica, 52(3), 701–720. doi:10.2307/1913472
• Santos Silva, J. M. C. and Tenreyro, S. (2006). The log of gravity. The Review of Economics and Statistics, 88(4), 641–658. doi:10.1162/rest.88.4.641