Count Regression Method Card Template (Statistics, Poisson and NB, R and Python)

Count Regression (Poisson and Negative Binomial): [Your study question or course]

A method card for regression models of counts: Poisson regression, rate ratios, offsets for unequal exposure time, the overdispersion check, quasi-Poisson, negative binomial regression (with the theta and alpha parameterisations used by R and Python) and a short note on zero inflation. It is for students preparing for a course or exam, thesis writers whose outcome is a number of events, and analysts who model incidents, visits, defects or complaints. It speeds up learning by putting the decisions in a fixed order (exposure, dispersion, zeros), working one data set fully in R and Python, then fading the support and testing you before you see the answers. The worked examples use fictional data; replace them with your own. Based on generalised linear model theory (Nelder and Wedderburn 1972), the NIST/SEMATECH e-Handbook of Statistical Methods (public domain) and the count-model references listed at the end.

Document control

When to use it

One-sentence rule: use count regression when the outcome is the number of times an event happened in a unit of time, space or exposure, and you want to know how the rate of events changes with the predictors.

Which count model?

The question and the data

Example question: at Harbour City Library (fictional), do branches with a children's area have a higher rate of damaged-item reports per opening hour than branches without one?

Design: observational, 5 branches x 4 months = 20 branch-months. The same branch appears four times; with only five branches we keep the example simple and treat rows as independent, and list the clustering as a limitation (see the mixed models card for the better model).

First look: the mean count is 9.8 and the variance is 51.4. Part of that variance comes from different opening hours and the predictor, so this is only a rough signal, but a variance five times the mean already suggests overdispersion.

The chart plots each branch-month's count against its opening hours. Longer opening goes with more reports, which is why the rate per hour, not the raw count, is the fair comparison.

Assumptions and how to check each

As a rule of thumb, a Pearson ratio up to about 1.2 to 1.5 is usually tolerated; above 2 the Poisson standard errors are clearly too small. The ratio is only reliable when the mean structure is right: a missing predictor or a wrong functional form also inflates it.

How it works

In plain language. A Poisson regression models the logarithm of the expected count as a straight line in the predictors. Exponentiated coefficients are rate ratios: the factor by which the expected number of events (per unit of exposure) is multiplied for a one-unit increase in a predictor. If rows are observed for different lengths of time, we add the log of the exposure to the model with its coefficient fixed at 1 (an offset), so the model compares rates per hour rather than raw counts. The Poisson distribution forces the variance to equal the mean; real counts usually vary more. The negative binomial model adds one parameter that lets the variance grow faster than the mean.

Poisson model with offset for row $i$ with exposure $t_i$:

so the rate per unit of exposure is

Mean and variance. Poisson: $\operatorname{E}(Y) = \operatorname{Var}(Y) = \mu$. Quasi-Poisson: $\operatorname{Var}(Y) = \phi\mu$, with $\phi$ estimated by the Pearson ratio. Negative binomial (NB2):

Dispersion statistic with $p$ coefficients:

R vs Python naming. MASS::glm.nb reports theta; statsmodels NegativeBinomial reports alpha. They are the same parameter inverted. In the example, theta = 5.53 in R and alpha = 0.1807 in Python, and 1 / 5.53 = 0.181. The statsmodels GLM family NegativeBinomial(alpha=...) treats alpha as fixed (default 1), not estimated, so use the discrete smf.negativebinomial (or estimate alpha first) unless you know alpha.

Zero inflation (note). Some data contain "structural" zeros from units that could never have the event (a branch closed for refurbishment reports no damage because it has no visitors). A zero-inflated Poisson or NB model mixes a logistic part (structural zero or not) with a count part; a hurdle model separates "any event" from "how many, given at least one". In R these are in the pscl package (zeroinfl(), hurdle()); in Python, statsmodels.discrete.count_model.ZeroInflatedPoisson and ZeroInflatedNegativeBinomialP. Fit them only when you have a reason for a separate zero process and the NB model still underpredicts zeros: an ordinary NB model often explains "too many zeros" on its own. Our library data contain no zeros at all, so zero inflation does not arise.

Worked example A (full): damaged items at Harbour City Library

Hypotheses. $H_0: \beta_{\text{kids}} = 0$ (rate ratio 1: the same rate of damaged-item reports per opening hour with and without a children's area) against $H_1: \beta_{\text{kids}} \neq 0$, two-sided, α = 0.05.

R code.

# install.packages("MASS")  # only if missing; MASS is a recommended package shipped with most R builds
library(MASS)

d <- data.frame(
  branch  = rep(c("Central", "Northgate", "Quayside", "Elm Park", "Riverside"), each = 4),
  hours   = c(300, 300, 280, 300, 180, 200, 180, 200, 220, 240,
              220, 240, 120, 140, 120, 140, 160, 160, 180, 180),
  kids    = rep(c(1, 0, 1, 0, 1), each = 4),
  damaged = c(22, 15, 30, 18, 4, 9, 2, 6, 7, 12,
              4, 6, 8, 3, 6, 12, 3, 10, 5, 14))

# 1. Poisson with exposure offset log(hours)
pois <- glm(damaged ~ kids + offset(log(hours)), family = poisson, data = d)
summary(pois)
exp(cbind(RR = coef(pois), confint(pois)))

# 2. Overdispersion check: Pearson chi-square / residual df
X2 <- sum(residuals(pois, type = "pearson")^2)
c(X2 = X2, df = df.residual(pois), ratio = X2 / df.residual(pois))

# 3. Quasi-Poisson: same estimates, SEs scaled by sqrt(dispersion)
qp <- glm(damaged ~ kids + offset(log(hours)), family = quasipoisson, data = d)
summary(qp)$coefficients

# 4. Negative binomial (theta parameterisation: Var = mu + mu^2 / theta)
nb <- glm.nb(damaged ~ kids + offset(log(hours)), data = d)
summary(nb)
exp(cbind(RR = coef(nb), confint(nb)))
c(theta = nb$theta, alpha = 1 / nb$theta)
AIC(pois, nb)

# Damage reports per 100 opening hours, from the NB model
exp(coef(nb)[1]) * 100                 # branches without a children's area
exp(sum(coef(nb))) * 100               # branches with a children's area

# The mistake: no offset
coef(glm(damaged ~ kids, family = poisson, data = d))

R output (R 4.5.1, MASS 7.3.65; output abridged: the section labels "Poisson:", "Quasi-Poisson:", "Negative binomial:" and "No offset:" and the bracketed notes after the two rates are added for reading; R does not print them).

Poisson:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.2426     0.1414 -22.929   <2e-16 ***
kids          0.2960     0.1639   1.806   0.0708 .  
    Null deviance: 57.157  on 19  degrees of freedom
Residual deviance: 53.741  on 18  degrees of freedom
AIC: 135.92
                   RR      2.5 %     97.5 %
(Intercept) 0.0390625 0.02921153 0.05091242
kids        1.3444604 0.98269311 1.87076198

       X2        df     ratio 
54.155864 18.000000  3.008659 

Quasi-Poisson:
              Estimate Std. Error    t value     Pr(>|t|)
(Intercept) -3.2425923  0.2452981 -13.218985 1.047193e-10
kids         0.2959928  0.2842160   1.041436 3.114602e-01

Negative binomial:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.2109     0.2059 -15.596   <2e-16 ***
kids          0.2330     0.2545   0.915     0.36    
              Theta:  5.53 
          Std. Err.:  2.84 
                    RR      2.5 %     97.5 %
(Intercept) 0.04032084 0.02685704 0.06073455
kids        1.26234287 0.76510662 2.08552660
    theta     alpha 
5.5326765 0.1807443 
     df      AIC
pois  2 135.9167
nb    3 123.4793

(Intercept) 
   4.032084            # per 100 hours, no children's area
[1] 5.089873           # per 100 hours, children's area

No offset:
(Intercept)        kids 
  1.8325815   0.6661185 

Python code.

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

d = pd.DataFrame({
    "hours":   [300, 300, 280, 300, 180, 200, 180, 200, 220, 240,
                220, 240, 120, 140, 120, 140, 160, 160, 180, 180],
    "kids":    [1] * 4 + [0] * 4 + [1] * 4 + [0] * 4 + [1] * 4,
    "damaged": [22, 15, 30, 18, 4, 9, 2, 6, 7, 12,
                4, 6, 8, 3, 6, 12, 3, 10, 5, 14]})
log_hours = np.log(d["hours"])

# 1. Poisson GLM with offset
pois = smf.glm("damaged ~ kids", data=d, family=sm.families.Poisson(),
               offset=log_hours).fit()
print(pois.summary().tables[1])
print("RR:", np.exp(pois.params).round(4).to_dict())

# 2. Overdispersion
print("Pearson chi2/df:", round(pois.pearson_chi2 / pois.df_resid, 3))

# 3. Quasi-Poisson: scale estimated from Pearson chi2 (z tests, not t as in R)
qp = smf.glm("damaged ~ kids", data=d, family=sm.families.Poisson(),
             offset=log_hours).fit(scale="X2")
print(qp.summary().tables[1])

# 4. Negative binomial (NB2) by maximum likelihood; alpha = 1 / theta
nb = smf.negativebinomial("damaged ~ kids", data=d, offset=log_hours).fit(disp=0)
print(nb.summary().tables[1])
print("alpha:", round(nb.params["alpha"], 4), " theta = 1/alpha:", round(1 / nb.params["alpha"], 3))
print("RR kids:", round(np.exp(nb.params["kids"]), 4),
      " Wald 95% CI:", np.exp(nb.conf_int().loc["kids"]).round(3).tolist())
print("AIC Poisson:", round(pois.aic, 3), " AIC NB:", round(nb.aic, 3))

Python output (Python 3.14, statsmodels 0.15.0; output abridged: the ==== table borders are omitted).

                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -3.2426      0.141    -22.929      0.000      -3.520      -2.965
kids           0.2960      0.164      1.806      0.071      -0.025       0.617
RR: {'Intercept': 0.0391, 'kids': 1.3445}
Pearson chi2/df: 3.009
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -3.2426      0.245    -13.219      0.000      -3.723      -2.762
kids           0.2960      0.284      1.041      0.298      -0.261       0.853
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -3.2109      0.208    -15.465      0.000      -3.618      -2.804
kids           0.2330      0.256      0.911      0.362      -0.268       0.734
alpha          0.1807      0.093      1.944      0.052      -0.001       0.363
alpha: 0.1807  theta = 1/alpha: 5.533
RR kids: 1.2623  Wald 95% CI: [0.765, 2.084]
AIC Poisson: 135.917  AIC NB: 123.479

Where R and Python differ, and why.

• The Poisson fits agree exactly. R's confint gives profile intervals for the rate ratio ([0.983, 1.871]); statsmodels gives Wald intervals on the log scale ([-0.025, 0.617], that is RR [0.975, 1.853]).
• Quasi-Poisson: same estimates and SEs, but R uses t tests on the residual df (p = .311) while statsmodels with scale="X2" uses z tests (p = .298).
• Negative binomial: the estimates and AIC agree; the SE of kids is 0.2545 in R and 0.256 in Python, because glm.nb takes the SEs from its final IRLS step (expected information, theta held at its estimate), while statsmodels uses the observed Hessian of the joint likelihood of the coefficients and alpha. R reports theta = 5.53, Python alpha = 0.1807 = 1/5.53. The NB CI in R is a profile interval ([0.765, 2.086]), in Python a Wald interval ([0.765, 2.084]).

Calculation by hand. With one binary predictor and an offset, the Poisson rate ratio equals the ratio of the two pooled rates:

• With a children's area: 146 reports in 2780 hours = 5.252 per 100 hours.
• Without: 50 reports in 1280 hours = 3.906 per 100 hours.
• Ratio $= 5.252 / 3.906 = 1.344 = e^{0.2960}$, matching the Poisson RR.
• Dispersion: $X^2 / df = 54.16 / 18 = 3.01$. The quasi-Poisson SE is the Poisson SE times $\sqrt{3.01} = 1.735$: $0.1639 \times 1.735 = 0.284$.

Decision. The Pearson ratio of 3.0 means the variance is about three times what the Poisson model allows, so the Poisson p = .071 and its CI are too optimistic. Both corrections widen the interval: quasi-Poisson p = .311, negative binomial p = .36. NB also has a much lower AIC (123.5 vs 135.9), so we report the NB model.

Effect size with CI. Rate ratio (NB) = 1.26, 95% CI [0.77, 2.09]: the estimated rate is 5.09 reports per 100 opening hours with a children's area and 4.03 without.

Interpretation. Branches with a children's area had an estimated 26% higher rate of damaged-item reports per opening hour, but the interval runs from 23% lower to 109% higher, so the data are compatible with no difference and with a doubling. The naive Poisson model nearly suggested a significant difference only because it ignored overdispersion. The analysis also treats four months from the same branch as independent; with five branches it cannot separate a "children's area" effect from other branch differences.

APA-style report. Damaged-item reports were modelled with negative binomial regression with log opening hours as an offset, because a Poisson model showed overdispersion (Pearson χ²/df = 3.01). The rate of reports was not clearly higher in branches with a children's area, rate ratio = 1.26, 95% CI [0.77, 2.09], z = 0.92, p = .360 (θ = 5.53).

Worked example B (half solved): rates from coefficients

The NB model gave $\hat{\beta}_0 = -3.2109$ and $\hat{\beta}_{\text{kids}} = 0.2330$, with log(hours) as the offset.

1. Rate per opening hour without a children's area: $e^{-3.2109} =$ [ ] reports per hour, which is [ ] per 100 hours.
2. Rate per opening hour with a children's area: $e^{-3.2109 + 0.2330} =$ [ ], which is [ ] per 100 hours.
3. Rate ratio: $e^{0.2330} =$ [ ]; check it equals the ratio of your answers to 2 and 1.
4. Expected reports in a month with 250 opening hours at a branch with a children's area: $250 \times$ [ ] $=$ [ ].
5. In R, theta = 5.53. The variance of the count when $\mu = 12.7$ is $\mu + \mu^2/\theta = 12.7 +$ [ ] $=$ [ ], compared with [ ] under Poisson.

Exercise C (unsolved)

Brightwell Bakery (fictional) counts customer complaints per week for 12 weeks. A new recipe was introduced in week 7. Weekly sales differ, so the exposure is thousands of loaves sold.

Tasks: (a) fit a Poisson model with offset(log(loaves_k)); (b) report the rate ratio for the new recipe with its 95% CI; (c) compute the Pearson dispersion ratio and decide whether NB is needed; (d) give the complaint rate per 1000 loaves before the change; (e) name one reason why this before-after comparison cannot prove the recipe caused the change.

Reading software output

Reporting template

Fill-in sentence. [Event] counts were modelled with [Poisson / negative binomial / quasi-Poisson] regression with log [exposure] as an offset ([reason, e.g. Pearson χ²/df = value]). Holding [other predictors] constant, [predictor] was associated with a rate ratio of [RR], 95% CI [lower, upper], z = [z], p [=/<] [p]; the estimated rate was [rate] per [unit of exposure] in [group] and [rate] in [reference group].

Note under the table: negative binomial, θ = [5.53]; offset log(opening hours); n = [20] branch-months.

Results checklist.

• Event and unit of observation defined
• Exposure and offset stated, or why none was needed
• Pearson dispersion ratio reported and the model choice justified
• Rate ratios with 95% CIs, CI type named
• At least one result as an absolute rate (per 100 hours, per 1000 loaves)
• Zeros checked against the model's expected zeros
• Clustering or repeated measures discussed
• Software, versions and the NB parameterisation (θ or α) named

Common mistakes and misconceptions

• Forgetting the offset. Without it, the example gives a "rate ratio" of $e^{0.666} = 1.95$ instead of 1.34, because branches with a children's area are also open longer. The model then compares counts, not rates.
• Putting log(exposure) in as an ordinary predictor when you mean a rate. That estimates a coefficient for it instead of fixing it at 1. (It can be a deliberate sensitivity check, but then say so.)
• Modelling the rate (count / hours) with linear regression. Rates from short exposures are noisy, the variance is not constant, and the model can predict negative rates.
• Skipping the dispersion check. The Poisson p = .071 looked "nearly significant"; after accounting for a dispersion of 3, p = .36.
• Confusing theta and alpha. A "dispersion of 5.53" in R and "0.18" in Python describe the same model. Large theta (small alpha) means close to Poisson.
• Using statsmodels GLM NegativeBinomial() without setting alpha. It fixes alpha at 1 instead of estimating it.
• Reaching for zero inflation because there are "many zeros". Many zeros are often just what a low mean or an NB model predicts. Compare observed and expected zeros first, and have a substantive reason for a separate zero process. The Vuong test for non-nested models is not a valid test of zero inflation against the non-inflated model, which is nested in it (Wilson 2015).
• Reading a rate ratio of 1.26 as "26 percentage points higher". It is a relative change in the rate: 26% higher.

Self-test questions

1. Why is the offset entered as log(hours) and not as hours?
2. A Poisson model gives Pearson χ² = 90 on 30 df. What is the dispersion ratio, by what factor are the quasi-Poisson SEs larger, and what would you do next?
3. R reports theta = 2 for a negative binomial model. What is alpha in statsmodels, and what is the variance of a count with mean 10?
4. A coefficient of -0.35 for weekend (1 = weekend day) in a model of library visits with an offset of log(opening hours). Interpret it as a rate ratio.
5. When would a zero-inflated model be justified for the number of books a resident borrowed last year?
6. Why can the overdispersion ratio be large even when the counts are Poisson given the right predictors?

My notes and next steps

Sources and further reading

Structure aligned with the method card used across this statistics library; formulas are standard generalised linear model results in our own notation.

• Nelder, J. A., and Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society, Series A, 135(3), 370-384.
• McCullagh, P., and Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall.
• Lawless, J. F. (1987). Negative binomial and mixed Poisson regression. Canadian Journal of Statistics, 15(3), 209-225.
• Lambert, D. (1992). Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics, 34(1), 1-14.
• Cameron, A. C., and Trivedi, P. K. (2013). Regression Analysis of Count Data (2nd ed.). Cambridge University Press.
• Venables, W. N., and Ripley, B. D. (2002). Modern Applied Statistics with S (4th ed.). Springer (the MASS package and glm.nb).
• Zeileis, A., Kleiber, C., and Jackman, S. (2008). Regression models for count data in R. Journal of Statistical Software, 27(8).
• Wilson, P. (2015). The misuse of the Vuong test for non-nested models to test for zero-inflation. Economics Letters, 127, 51-53.
• NIST/SEMATECH e-Handbook of Statistical Methods, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
• Further reading by name: de Smith, M. J., Statistical Analysis Handbook (StatsRef, 2024 edition); OpenIntro Statistics; APA Publication Manual (7th ed.) and APA JARS-Quant.

Answer key

Example B.

1. $e^{-3.2109} = 0.04032$ per hour, 4.03 per 100 hours.
2. $e^{-2.9779} = 0.05090$ per hour, 5.09 per 100 hours.
3. $e^{0.2330} = 1.262$; $5.090 / 4.032 = 1.262$.
4. $250 \times 0.05090 = 12.7$ reports.
5. $12.7^2 / 5.53 = 29.2$, so variance $= 12.7 + 29.2 = 41.9$, compared with 12.7 under Poisson.

Exercise C (R 4.5.1; output abridged and relabelled: the values come from coef(summary()), exp(cbind(coef(), confint())) and the Pearson ratio, with the variable named new_recipe):

              Estimate Std. Error   z value     Pr(>|z|)
(Intercept)  0.8788643  0.1301889  6.750685 1.471488e-11
new_recipe  -0.4348612  0.2033209 -2.138792 3.245250e-02
rate ratio and profile 95% CI:
(Intercept) 2.4081633 [1.8447918, 3.0759689]
new_recipe  0.6473545 [0.4319296, 0.9609721]
Pearson chi2 / df: 0.1718832

(b) Rate ratio = 0.65, 95% CI [0.43, 0.96], z = -2.14, p = .032: about 35% fewer complaints per 1000 loaves after the change. (c) The ratio is 0.17, well below 1 (underdispersion: the counts are more regular than Poisson predicts). NB is not needed: it can only add variance, never remove it. If the counts really are underdispersed, the Poisson SEs are too large, so the Poisson test is conservative. With only 12 weeks the ratio is itself very imprecise, so do not shrink the SEs with a quasi-Poisson scale below 1; report the Poisson model and note the low dispersion. (d) $e^{0.8789} = 2.41$ complaints per 1000 loaves. (e) Anything else that changed over time (season, staff, customer mix, a new complaint form) is confounded with the recipe; there is no concurrent control group.

Self-test.

1. The model is for the log of the expected count. Since $\log(\text{rate} \times t) = \log(\text{rate}) + \log t$, adding $\log t$ with a coefficient of 1 turns the model for counts into a model for rates per unit of exposure.
2. $90 / 30 = 3.0$; SEs are $\sqrt{3} = 1.73$ times larger. Check the mean model (missing predictors, non-linearity), then use quasi-Poisson or a negative binomial model.
3. $\alpha = 1/2 = 0.5$; variance $= 10 + 0.5 \times 100 = 60$.
4. $e^{-0.35} = 0.70$: on weekend days the visit rate per opening hour is about 30% lower than on weekdays, other predictors held constant.
5. When some residents could never borrow (for example, people who are not library members) and you cannot identify them in the data, so a separate "structural zero" process is plausible, and the NB model still predicts fewer zeros than observed.
6. A wrong mean model (a missing important predictor, a missing interaction, wrong functional form) or outliers make the Pearson residuals large, which inflates the ratio even when the true conditional distribution is Poisson.

Revision history