← All examples

Count Models

Regression for count data, threaded by one problem — extra-Poisson variation, the near-universal fact that real counts vary more than a plain Poisson allows, and the different ways to model it. It starts with the two core GLMs, Poisson and Negative Binomial regression, where the NB's dispersion parameter restores calibrated inference under overdispersion. The hierarchical Poisson GLMM adds group random effects on the log-rate, and the hierarchical binomial GLMM is its twin on grouped proportions — together completing the Gaussian / Poisson / binomial hierarchical-GLM grid. Zero-inflated Poisson and Negative Binomial handle counts with far more zeros than any Poisson can produce, splitting structural from sampling zeros via a data-augmentation Gibbs sampler. And the Besag–York–Mollié spatial model closes the set — disease mapping where an intrinsic-CAR prior pools each area toward its map neighbours rather than a global mean. Each is built from scratch and cross-checked across PyMC and R.

What makes a count different. A count is a whole number of things that happened — seizures in a fortnight, seeds that germinated on a plate, cases of a disease in a district. Three properties separate it from an ordinary measurement. It cannot be negative. It has no natural ceiling (or, for grouped counts, a known one: you cannot germinate more seeds than you planted). And — the property that drives this entire section — its variability is tied to its level rather than free of it. A district averaging 2 cases will never show a deviation of ±50; one averaging 500 will never show a deviation of ±2. Ordinary regression assumes the spread of the errors is the same everywhere, and for counts that assumption is not merely wrong, it is impossible. The Poisson builds the tie in exactly: variance = mean. That is equidispersion, it is the strongest assumption made anywhere in this section, and every model below is a response to the fact that real counts vary more than it permits.

GLM, and what the link does. Every model here is a generalised linear model — the extension of ordinary regression to outcomes that are not Normal. The machinery is the same in each case: form a linear combination of the predictors, then pass it through a link function that maps it onto the scale the outcome actually lives on. Counts cannot be negative, so the link is the log: the model predicts log(rate), and exponentiating returns a rate that is positive whatever the coefficients do. When a random effect per group is added on top, the model becomes a GLMM, the second M standing for mixed — mixing coefficients that are fixed across the whole sample with coefficients that vary from group to group. Four of the five examples below are GLMMs.

Why the log link makes everything a ratio. A consequence worth stating, because it governs how every result on these pages is reported. On a log scale, adding a coefficient multiplies the rate: a coefficient of −0.23 means the rate is multiplied by e^−0.23 = 0.79, a 21% reduction. That multiplier is the rate ratio, and it is what count models report instead of a difference, because a difference would mean something different for a hospital doing 47 operations than for one doing 810. The same object appears as a relative risk in disease mapping and as an incidence rate ratio in epidemiology; they are the same number under different names.

The offset, which is how counts become rates. Comparing raw counts across units of different size is meaningless — a district with 39 lip-cancer cases and one with 9 are not comparable until you know how many people are in each. The standard device is an offset: a term added to the linear predictor whose coefficient is fixed at 1 rather than estimated. Putting log(expected cases) in as an offset is exactly equivalent to modelling the ratio of observed to expected, so the coefficients become statements about rates while the likelihood stays a proper count likelihood. In the spatial model below the offset is age-standardised expected cases; elsewhere it would be population, exposure time, or number of trials.

One construction, four kinds of outcome. The hierarchical models in this collection are worth seeing as a single grid rather than as separate examples, because they are the same construction with the likelihood swapped. In every case each unit — a rat, a patient, a plate, a district — gets its own intercept, those intercepts are treated as draws from one common Normal distribution, and that distribution’s spread is estimated alongside everything else. What changes between them is only what sits on top. With a Gaussian outcome it is the Rats growth model, in the regression section. With a Poisson outcome it is Epil, below. With a binomial outcome it is Seeds, also below. With a binary outcome measured repeatedly it is the toenail trial, in limited dependent variables. Same random intercept, four likelihoods — and the only thing the choice of likelihood really changes is the sampler: conjugate and exact for the Gaussian, Pólya–Gamma for the binary and binomial logits, and Metropolis-within-Gibbs for the Poisson, which has no such trick available.