Poisson regression for counts

Setup#

pip install scikit-learn

Outputs verified with scikit-learn 1.7.2 on CPU.

Adverts and shop visits#

import numpy as np
from sklearn.linear_model import LinearRegression, PoissonRegressor

rng = np.random.default_rng(3)
ads = rng.integers(0, 6, 60).astype(float)          # adverts run each week
visits = rng.poisson(np.exp(0.8 + 0.35 * ads))      # customer visits that week

X = ads.reshape(-1, 1)
lin = LinearRegression().fit(X, visits)
poi = PoissonRegressor(alpha=0).fit(X, visits)

print(f"linear : visits = {lin.intercept_:.1f} + {lin.coef_[0]:.1f} * ads")
print(f"poisson: each advert multiplies visits by {np.exp(poi.coef_[0]):.2f}")
for n in (0, 5, 10):
    a = np.array([[float(n)]])
    print(f"ads={n:2d}: linear {lin.predict(a)[0]:6.1f}   poisson {poi.predict(a)[0]:6.1f}")

linear : visits = 0.8 + 2.5 * ads
poisson: each advert multiplies visits by 1.49
ads= 0: linear    0.8   poisson    2.1
ads= 5: linear   13.2   poisson   15.2
ads=10: linear   25.7   poisson  111.8

The walkthrough#

The models tell different stories. Linear: "each advert adds 2.5 visits, whatever the current level." Poisson: "each advert multiplies visits by 1.49." The data was generated multiplicatively (each advert multiplies the rate by exp(0.35) ≈ 1.42), and only one model can even say that sentence.

At zero adverts, linear predicts 0.8 visits — beneath the true quiet-week average of about 2.2, dragged low by forcing one straight line through curved data. Add a feature with a negative effect (rain, roadworks) and linear regression will cross into negative visits without blinking.

At ten adverts, the models disagree by 4x (25.7 vs 111.8). Both are extrapolating beyond the training range of 0–5, so treat both numbers with suspicion — but note the shapes: exponential growth versus a straight line is a modelling decision with enormous consequences outside the data. Exponential extrapolation can also overshoot badly; honesty about the range matters in both directions.

np.exp(poi.coef_[0]) is the reading ritual. Poisson coefficients live in log-space; exponentiate to get multipliers. A coefficient of -0.5 means exp(-0.5) ≈ 0.61 — that feature cuts the rate by 39%.

alpha=0 switches off PoissonRegressor's default ridge penalty (alpha=1.0) so the comparison with plain OLS is fair. In production, keep some regularisation.

Common mistakes#

Negative or fractional targets. PoissonRegressor refuses negative targets outright — ValueError: Some value(s) of y are out of the valid range of the loss 'HalfPoissonLoss'. Rates like "visits per opening hour" (fractional) belong in the exposure pattern below, not in the target.

Ignoring exposure. Comparing weekly counts across shops open different hours is comparing apples to weeks. Model counts with the time-open as an offset — in statsmodels, sm.GLM(..., family=sm.families.Poisson(), offset=np.log(hours_open)) — so coefficients describe rates, not raw totals.

Trusting Poisson's variance promise. The Poisson family assumes wobble equals the mean: an average of 5 implies variance 5. Real counts are usually messier — one viral post, one festival — a state called overdispersion. Your rate estimates stay reasonable, but the error bars become too confident. Check, and escalate to negative binomial if needed (researcher block).

Log-transforming the target instead. log(visits + 1) into OLS looks similar and behaves worse: the +1 is arbitrary, zeros distort, and predictions are biased on the way back. Poisson regression is the correct version of that instinct.

Try it yourself#

Generate a second dataset with visits = rng.poisson(np.exp(0.8 + 0.35 * ads)) * rng.integers(1, 4, 60) — bursty, overdispersed arrivals — and refit. Compare the mean and variance of visits at each ads level to see the broken promise. Do the rate estimates survive?

What to learn next#

• Generalised linear models — the machine this attachment clicks into.
• Quantile regression — when you need the busy-day tail, not the average.
• Time series basics — counts over time bring their own structure.