# 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 | Field | Value | |---|---| | Course or project | [STAT 310 Generalised linear models / library service evaluation] | | Learner or analyst | [Your name] | | Data set | [*Example:* Harbour City Library damaged-item reports, 5 branches x 4 months, fictional] | | Software and versions | [*Example:* R 4.5.1 with MASS 7.3.65; Python 3.14 with numpy 2.5.4, pandas 3.1.0, statsmodels 0.15.0] | | Date | YYYY-MM-DD | | Version | [0.1] | ## 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. | Use when | Do not use when | |---|---| | The outcome is a non-negative integer count (incidents, visits, defects) | The outcome is "k successes out of n trials" with known n: use binomial (logistic) regression | | Rows may have different exposure (opening hours, months, population): use an offset | The counts are so large and far from zero that a linear model on a log or square root scale behaves well, and you prefer it for simplicity | | You want rate ratios adjusted for other predictors | Counts are repeated over time for the same units and you care about the time pattern: consider mixed models or time series models | | Counts are independent given the predictors | Rows are clustered (same branch in several months) and the clustering is strong: add random effects or cluster-robust SEs | **Which count model?** | Situation | Model | |---|---| | Variance about equal to the mean (Pearson χ²/df near 1) | Poisson | | Variance larger than the mean, you only need valid SEs | Quasi-Poisson, or Poisson with robust (sandwich) SEs | | Variance larger than the mean, you want a full likelihood (AIC, predictions of counts) | Negative binomial (NB2) | | More zeros than the Poisson or NB model predicts, from a separate "never" process | Zero-inflated or hurdle 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? | Variable | Role | Measurement level | Unit or coding | |---|---|---|---| | `damaged` | Outcome | Count | Damaged-item reports in the month | | `hours` | Exposure (offset) | Continuous (ratio) | Opening hours in the month; enters as $\log(\text{hours})$ | | `kids` | Predictor | Binary | 1 = branch has a children's area, 0 = no (reference) | | `branch` | Identifier | Nominal | Five branches, four months each | **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). | Row | Branch | hours | kids | damaged | Row | Branch | hours | kids | damaged | |---|---|---|---|---|---|---|---|---|---| | 1 | Central | 300 | 1 | 22 | 11 | Quayside | 220 | 1 | 4 | | 2 | Central | 300 | 1 | 15 | 12 | Quayside | 240 | 1 | 6 | | 3 | Central | 280 | 1 | 30 | 13 | Elm Park | 120 | 0 | 8 | | 4 | Central | 300 | 1 | 18 | 14 | Elm Park | 140 | 0 | 3 | | 5 | Northgate | 180 | 0 | 4 | 15 | Elm Park | 120 | 0 | 6 | | 6 | Northgate | 200 | 0 | 9 | 16 | Elm Park | 140 | 0 | 12 | | 7 | Northgate | 180 | 0 | 2 | 17 | Riverside | 160 | 1 | 3 | | 8 | Northgate | 200 | 0 | 6 | 18 | Riverside | 160 | 1 | 10 | | 9 | Quayside | 220 | 1 | 7 | 19 | Riverside | 180 | 1 | 5 | | 10 | Quayside | 240 | 1 | 12 | 20 | Riverside | 180 | 1 | 14 | 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. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Damaged-item reports against opening hours by children's area", "width": 380, "height": 240, "data": {"values": [ {"hours": 300, "damaged": 22, "area": "Children's area"}, {"hours": 300, "damaged": 15, "area": "Children's area"}, {"hours": 280, "damaged": 30, "area": "Children's area"}, {"hours": 300, "damaged": 18, "area": "Children's area"}, {"hours": 180, "damaged": 4, "area": "No children's area"}, {"hours": 200, "damaged": 9, "area": "No children's area"}, {"hours": 180, "damaged": 2, "area": "No children's area"}, {"hours": 200, "damaged": 6, "area": "No children's area"}, {"hours": 220, "damaged": 7, "area": "Children's area"}, {"hours": 240, "damaged": 12, "area": "Children's area"}, {"hours": 220, "damaged": 4, "area": "Children's area"}, {"hours": 240, "damaged": 6, "area": "Children's area"}, {"hours": 120, "damaged": 8, "area": "No children's area"}, {"hours": 140, "damaged": 3, "area": "No children's area"}, {"hours": 120, "damaged": 6, "area": "No children's area"}, {"hours": 140, "damaged": 12, "area": "No children's area"}, {"hours": 160, "damaged": 3, "area": "Children's area"}, {"hours": 160, "damaged": 10, "area": "Children's area"}, {"hours": 180, "damaged": 5, "area": "Children's area"}, {"hours": 180, "damaged": 14, "area": "Children's area"} ]}, "mark": {"type": "point", "filled": true, "size": 70}, "encoding": { "x": {"field": "hours", "type": "quantitative", "title": "Opening hours in the month", "scale": {"zero": false}}, "y": {"field": "damaged", "type": "quantitative", "title": "Damaged-item reports"}, "color": {"field": "area", "type": "nominal", "title": "Branch"} } } ``` ## Assumptions and how to check each | Assumption | How to check | R | Python | If violated | |---|---|---|---|---| | Counts of independent events | Study design: no contagion between events, no clustering ignored | Check the design | Check the design | Random effects (`lme4::glmer(..., family = poisson)`), cluster-robust SEs, or NB | | Log-linear mean: log of the rate is linear in the predictors | Add a squared or spline term; plot Pearson residuals against each predictor | `anova(m, update(m, . ~ . + I(x^2)), test = "LRT")` | Compare `llf` of nested fits | Transform the predictor, add a spline | | Exposure handled correctly | Each row's exposure known and entered as `offset(log(exposure))` | `offset(log(hours))` in the formula | `offset=np.log(d["hours"])` (or `exposure=d["hours"]`) | Add the offset; never divide the count by exposure and model the ratio as a count | | Equidispersion (variance = mean) | Pearson χ² / residual df near 1; well above 1 is overdispersion | `sum(residuals(m, "pearson")^2) / df.residual(m)` | `m.pearson_chi2 / m.df_resid` | Quasi-Poisson, robust SEs or negative binomial | | Zeros consistent with the model | Compare observed zeros with $\sum_i P(Y_i = 0)$ from the fitted model | `sum(dpois(0, fitted(m)))` vs `sum(d$y == 0)` | `stats.poisson.pmf(0, m.mu).sum()` | NB often fixes excess zeros; otherwise a zero-inflated or hurdle model | | No overly influential rows | Cook's distance for GLMs | `plot(m, which = 5)` | `m.get_influence().cooks_distance` | Check data, report with and without | 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$: $$ Y_i \sim \text{Poisson}(\mu_i), \qquad \log \mu_i = \log t_i + \beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik} $$ so the rate per unit of exposure is $$ \frac{\mu_i}{t_i} = \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik}), \qquad \text{RR}_j = e^{\beta_j} $$ **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): $$ \operatorname{Var}(Y) = \mu + \frac{\mu^2}{\theta} = \mu + \alpha\mu^2, \qquad \alpha = \frac{1}{\theta} $$ **Dispersion statistic** with $p$ coefficients: $$ \hat{\phi} = \frac{X^2}{n - p} = \frac{1}{n - p}\sum_{i=1}^{n} \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i} $$ | Symbol | Meaning | |---|---| | $Y_i$, $\mu_i$ | Observed count and its expected value for row $i$ | | $t_i$ | Exposure (opening hours); $\log t_i$ is the offset | | $\beta_0$ | Log of the event rate per unit of exposure at all predictors 0 | | $e^{\beta_j}$ | Rate ratio for a one-unit increase in $x_j$, others held constant | | $X^2$ | Pearson chi-square statistic | | $\hat{\phi}$ | Estimated dispersion (1 for a true Poisson) | | $\theta$ | NB "size" parameter in R (`glm.nb`); large $\theta$ means close to Poisson | | $\alpha$ | NB dispersion parameter in statsmodels; $\alpha = 1/\theta$, and $\alpha = 0$ is Poisson | **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. ```mermaid flowchart TD A["Outcome is a count per unit"] --> B{"Exposure differs between rows?"} B -->|Yes| C["Add offset log of exposure"] B -->|No| D["Fit Poisson"] C --> D D --> E["Pearson chi-square divided by df"] E -->|"About 1"| F["Keep Poisson"] E -->|"Well above 1"| G["Check the mean model first"] G --> H{"Need a full likelihood?"} H -->|No| I["Quasi-Poisson or robust SEs"] H -->|Yes| J["Negative binomial"] F --> K{"Excess zeros with a reason?"} I --> K J --> K K -->|Yes| L["Zero-inflated or hurdle model"] K -->|No| M["Report rate ratios with CIs"] L --> M ``` ## 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.** ```r # 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). ```text 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.** ```python 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). ```text 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. | Week | loaves_k | new_recipe | complaints | Week | loaves_k | new_recipe | complaints | |---|---|---|---|---|---|---|---| | 1 | 4.0 | 0 | 9 | 7 | 4.4 | 1 | 6 | | 2 | 4.2 | 0 | 11 | 8 | 4.6 | 1 | 8 | | 3 | 3.8 | 0 | 8 | 9 | 4.0 | 1 | 5 | | 4 | 4.5 | 0 | 12 | 10 | 4.3 | 1 | 9 | | 5 | 4.1 | 0 | 10 | 11 | 4.8 | 1 | 7 | | 6 | 3.9 | 0 | 9 | 12 | 4.2 | 1 | 6 | 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 | Output line | Where | Meaning | Report it? | |---|---|---|---| | `Estimate` / `coef` | R and Python | Change in log rate per unit; exponentiate for the rate ratio | As RR with CI | | `(Intercept)` | R and Python | Log rate per unit of exposure at all predictors 0 | As a baseline rate if meaningful | | `Residual deviance ... on 18 degrees of freedom` | R | Deviance-based fit; deviance/df is a second dispersion check | Not usually | | Pearson χ² / df | computed | Dispersion; about 1 for Poisson | Yes, to justify the model | | `Dispersion parameter for quasipoisson family` | R `summary(qp)` | $\hat{\phi}$, here 3.01 | Yes, with quasi-Poisson | | `Theta` and `Std. Err.` | R `glm.nb` | NB size parameter; small theta = strong overdispersion | Yes, θ | | `alpha` | Python `negativebinomial` | NB dispersion, $1/\theta$; a test of alpha = 0 is a boundary test, so its p-value is conservative | Yes, α (or convert to θ) | | `AIC` | R and Python | Compare Poisson and NB on the same data (same package) | When choosing the model | ## 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]. | Predictor | *b* | SE | Rate ratio | 95% CI | *p* | |---|---|---|---|---|---| | *Example:* Intercept (rate per hour) | -3.211 | 0.206 | 0.040 | [0.027, 0.061] | < .001 | | *Example:* Children's area (vs none) | 0.233 | 0.255 | 1.26 | [0.77, 2.09] | .360 | | [Predictor] | [ ] | [ ] | [ ] | [ ] | [ ] | 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 | Item | Entry | |---|---| | Confidence with this method (1-5) | [3] | | What I can now do | [Fit Poisson and NB with an offset and justify the choice] | | Still unclear | [Hurdle vs zero-inflated models] | | Mistakes to log in academic/stats-error-log | [Forgot the offset] | | Next review date (spacing) | YYYY-MM-DD | | Next template | [academic/stats-mixed-models] | ## 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`): ```text 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 | Version | Date | Author | Change | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Name] | First draft from the template | | [0.2] | YYYY-MM-DD | [Name] | [Switched from Poisson to negative binomial] |