A Poisson regression is a generalized linear model with a log link, used when your outcome is a non-negative integer count of events, such as complaints per month, hospital visits per year, or accidents per period of service. Running one takes four moves: describe the count, fit the model with a log link, test whether the variance really equals the mean, then exponentiate the coefficients to read incidence rate ratios. It takes maybe fifteen minutes in R or SPSS once the data are tidy, and the rest of this guide is the workflow, not the theory.
I use one small worked example throughout so you can compare output across tools: complaints recorded at 24 retail branches in a single month. The count is complaints, the numeric predictor is staff, and the categorical predictor is region with North as the reference group. Every number below is illustrative, so treat the figures as teaching examples rather than findings.
Table of Contents
- 1What You Need
- 2Step-by-Step: How to Run Poisson Regression for Count Data
- 31. Define the count outcome and predictors
- 42. Inspect and prepare the data
- 53. Choose and fit a Poisson regression model
- 64. Check the assumptions and model fit
- 75. Interpret the coefficients and predicted counts
- 86. Report the Poisson regression results
- 9Common Mistakes
- 10Frequently Asked Questions
- 11What data can I use with Poisson regression?
- 12How do I test for overdispersion in a Poisson model?
- 13When should I use negative binomial instead of Poisson regression?
- 14How do I handle zero-inflated count data?
- 15Which software should I use, and how do I read the output?
- 16Conclusion
What You Need
Before opening any software, confirm that five things are true of your dataset.
- The outcome is a genuine count. Values are 0, 1, 2, 3 with no decimals, and each unit is one whole event. A percentage, a rating on a 1-to-10 scale, or a per-person average is not a count, and forcing a count model onto it produces numbers nobody can interpret.
- The predictors are the ones you actually need. Numeric predictors enter as continuous values; categorical predictors enter as factors with one level set as the reference group.
- You know the unit of analysis. One row is one branch per month, one patient, one trip. If a patient contributes three visits, that is either three rows or a clustered model, not one row with a count of three and no attention to the repeat structure.
- You know whether exposure varies. If branches were open for different numbers of days, or people were followed for different lengths of time, you need an offset. Without it you are comparing raw totals and attributing the difference to your predictors.
- You have software. R with MASS and pscl, Python with statsmodels, SPSS with the Generalized Linear Models menu, or Stata. All four fit the same model and produce comparable numbers.
Write the research question in one sentence first. For my example: does the number of complaints per branch-month change with staffing level and region?
Step-by-Step: How to Run Poisson Regression for Count Data
1. Define the count outcome and predictors
Start by separating your outcome into two parts: the count itself, and the opportunity for events to happen. In the complaints data, the count is complaints per branch-month and the exposure is days open.
Most of the questions people get stuck on show up here, not in the software. A researcher whose outcome is 2.4 visits per month is looking at a rate, not a count. Two options work: model the underlying count with an offset for exposure, or model the rate directly with a continuous model such as a log-normal or gamma regression. Rounding 2.4 to 2 throws away information and can bias the results.
Then check that the outcome actually behaves like a count. A quick histogram answers most of it.
hist(branches$complaints)
mean(branches$complaints) # expected count
var(branches$complaints) # observed variance
table(branches$complaints == 0) # how many zeros
If the mean and variance are wildly different, stop and read step 4 before fitting anything.
2. Inspect and prepare the data
Count models drop every row with a missing value, silently. With 24 branches that is invisible; with 3,000 rows it can quietly shrink your sample. Check first, then decide.
colSums(is.na(branches)) # missing per column
branches <- branches[complete.cases(branches), ] # or model them explicitly
duplicated(branches) # repeated branch-month rows
Other preparation worth doing now rather than later:
- Drop duplicate records. Two identical branch-month rows usually mean a merge went wrong.
- Check that the count has no negatives. A value of -1 means the variable is a difference, not a count.
- Set the reference level of every categorical predictor deliberately, and write down which level it is. You need it for the results table later.
- Look at extreme counts. A single branch with 60 complaints against a mean of 4 will swing the fit; check whether it is a data error or a genuinely different case.
3. Choose and fit a Poisson regression model
The model fits the log of the expected count: log(mu) = b0 + b1*X1 + … + bk*Xk. Parameters are estimated by maximum likelihood. The log link is what makes the exponentiated coefficients multiplicative instead of additive, and that is the entire reason for choosing this family over Gaussian regression.
In R, glm() with family = poisson does it:
library(MASS)
m1 <- glm(complaints ~ staff + region, data = branches, family = poisson)
summary(m1)
exp(cbind(IRR = coef(m1), conf.int(m1)))
In Python, statsmodels takes a formula and a family object:
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
m1 = smf.glm("complaints ~ staff + C(region)", data=branches,
family=sm.families.Poisson()).fit()
print(m1.summary())
print(np.exp(m1.params))
If you want penalised estimates for a wide predictor set, PoissonRegressor with an alpha value takes the same formula and adds ridge or lasso shrinkage. Stick with the plain GLM while you are learning.
In SPSS, the path is Analyze > Generalized Linear Models > Poisson. Put complaints in Dependent, region in Factors, staff in Covariates. That is the whole dialog. Under Options, the default fixed Scale of 1 gives a standard Poisson fit; switching it to Pearson Chi-Square gives you the quasi-Poisson dispersion and sandwich standard errors instead.
In Stata, one line:
poisson complaints staff i.region
estimates, store(P)
nbreg complaints staff i.region
estimates table P ., stats(N) b(%4.2f) se(%4.2f) p(%4.2f)
For rate data, add the exposure on the right-hand side in all four tools. In R that is offset(log(opdays)), in SPSS it is a value entered on the Offset tab, and in Stata it is the exposure() option:
poisson complaints staff i.region, exposure(opdays)
Interactions work as you would expect, using : in R and Python, and ## or c.staff#c.staff in Stata.
4. Check the assumptions and model fit
The Poisson model assumes the variance of the count equals its mean. That is the assumption real count data breaks most often, and a reader who skips this step will report standard errors that are too small.
Overdispersion. Compare mean and variance, then compute the dispersion directly from the fitted model:
phi <- sum(residuals(m1, type = "pearson")^2) / df.residual(m1)
phi # values well above 1 mean variance is too large for the Poisson model
In Stata, estat gof prints the same statistic with a p-value. In SPSS, set Scale to Pearson Chi-Square and read the dispersion line in the Iterations table. A dispersion near 1 is what the Poisson model expects; anything near 2 or higher is overdispersion.
Residual structure. Plot Pearson or deviance residuals against fitted values. A clean Poisson fit shows an unstructured band centred on zero. Fan shapes mean variance depends on the mean. A smooth curve means the relationship is non-linear on the log scale and a spline or a quadratic term belongs in the model.
plot(m1, which = 1) # residuals vs fitted
qqplot(residuals(m1, type = "pearson"))
Zero inflation. Compare the observed share of zeros with what the fitted model predicts. In the example, if 9 of 24 branches recorded zero complaints while the fitted model expects 4, that is 21 percentage points more zeros than the model accounts for, enough to distort a mean.
Independence. Plot counts against time order, by site, or by subject. Repeated measures per patient, per school or per site break independence and make standard errors too small. The fixes are cluster-robust standard errors, a mixed model, or GEE with a correlation structure.
Here is the decision rule I use once the diagnostics are in hand:
| Model | Use when | Key assumption | R | Stata |
|---|---|---|---|---|
| Poisson | Variance roughly equals the mean, no excess zeros | Equidispersion | glm(family = poisson) | poisson |
| Quasi-Poisson | Mild overdispersion, you want the Poisson coefficients with corrected standard errors | Mean model correct, variance free | glm(family = quasipoisson) | poisson, vce(robust) |
| Negative binomial | Substantial overdispersion, or variance grows faster than the mean | Variance = mu + mu squared times kappa | glm.nb() | nbreg |
| Zero-inflated Poisson | Too many zeros plus overdispersion from two generating processes | Mixture of a degenerate zero and a Poisson | pscl::zeroinfl | zip |
| Hurdle | Too many zeros where the zeros themselves are meaningful | Two-part: any event, then size | pscl::hurdle | crreg (counts) |
A likelihood ratio test between the Poisson and negative binomial fits settles the choice when the evidence is close. In R that is anova(m1, m2, test = "Chisq"). If the two models disagree about significance, that disagreement is a finding worth reporting, not something to hide.
5. Interpret the coefficients and predicted counts
Exponentiate. The Estimate column on the log scale becomes an incidence rate ratio, the multiplicative change in the expected count for a one-unit increase in the predictor, holding other variables fixed.
| Output column | What it means |
|---|---|
| Estimate (B) | Change in the log of the expected count. Not readable on its own. |
| Std. Error | Standard error of B on the log scale. |
| Wald z | Estimate divided by standard error. Values beyond plus or minus 1.96 are significant at p < .05. |
| p-value | Two-sided test that the log rate equals zero. |
| IRR | e raised to B. The rate multiplier. Above 1 means a higher expected count. |
| 95% CI | Confidence limits for the IRR. If it crosses 1, the effect is not distinguishable from no change. |
Worked reading for the example: if staff has B = 0.024, then IRR = e^0.024 = 1.024. One additional staff member is associated with about a 2.4 percent higher expected number of complaints. On the response scale that is a rate difference, not a probability change, and saying anything about probability would be wrong.
For region, North is the reference, so an IRR of 1.34 for West means the fitted expected count in West branches is 34 percent higher than in North branches, on the multiplicative scale. Whether that difference is trustworthy depends entirely on the confidence interval, as step 6 shows. The reference category itself never appears in the table, which is why you write it down in step 2.
The intercept is rarely of interest. It is the expected count when every predictor sits at zero, and a staff count of zero usually makes that prediction meaningless. Centering your predictors at a plausible value gives you a usable intercept:
m1 <- glm(complaints ~ staff_c + region, data = branches,
family = poisson, offset(log(branches$opdays)))
predict(m1, type = "response") # counts, back on the original scale
For a reader who wants probabilities, convert the fitted mean into a probability distribution with dpois(). That step is about presentation, not about the model.
6. Report the Poisson regression results
Report four things: the model and its link, the variables and reference categories, the fit evidence including the dispersion statistic, and the rate ratios with confidence intervals. An APA 7 table carries five columns.
| Predictor | B | SE | Wald z | p | IRR | 95% CI |
|---|---|---|---|---|---|---|
| Intercept | -1.86 | 0.31 | -6.00 | <.001 | 0.16 | [0.09, 0.29] |
| Staff (per person) | 0.024 | 0.009 | 2.67 | .008 | 1.02 | [1.01, 1.04] |
| Region: South | 0.11 | 0.24 | 0.46 | .64 | 1.12 | [0.70, 1.79] |
| Region: West | 0.29 | 0.23 | 1.26 | .21 | 1.34 | [0.85, 2.10] |
A template results paragraph you can adapt:
A Poisson regression was fitted to predict the number of complaints per branch-month from staffing level and region (North as reference), with log days open included as an offset. The model was estimated by maximum likelihood. A dispersion parameter of 1.08 indicated no substantial overdispersion, so the Poisson specification was retained. Each additional staff member was associated with a 2.4 percent increase in the expected number of complaints, IRR = 1.02, 95% CI [1.01, 1.04], z = 2.67, p = .008. The point estimate for West branches was a 34 percent higher rate than North branches, IRR = 1.34, 95% CI [0.85, 2.10], but the interval crossed 1, so that difference was not statistically distinguishable from no change. Results are illustrative.
Note the last clause. Always say whether the model is Poisson, quasi-Poisson or negative binomial, and always state the dispersion or the test that settled it.
Common Mistakes
Running ordinary least squares on counts. Predictions can go negative, errors are not normal, and variance does not grow with the mean. Use a Poisson GLM, or a negative binomial if the variance is much larger.
Forgetting the log link when reporting. A coefficient of 0.024 is not a 0.024 increase in complaints. Exponentiate it, and describe the result as a percentage change in the expected count.
Ignoring overdispersion. This is the single most common flaw in published count analyses. Compute the dispersion, and switch to quasi-Poisson or negative binomial when it runs well above 1.
Adding the exposure as an ordinary covariate. An offset has its coefficient fixed at 1. Enter it as a covariate and the software will estimate a coefficient for it, which changes the model you meant to fit.
Using percentages and rates as if they were counts. If the outcome has decimals, it is a rate or a proportion. Model the underlying count with an exposure offset, or choose a model for the continuous outcome.
Leaving repeated observations unlabelled. Several rows per student, patient or site violate independence. Use cluster-robust standard errors, GEE, or a mixed model.
Reading the intercept as an average. It is the expected count at zero on every predictor. If zero is outside the range of your data, center the predictors or ignore the intercept.
Skipping the zero check. Excess zeros shift the mean and the standard errors. Compare the observed zeros with the fitted zeros before deciding whether a zero-inflated or hurdle model fits.
Two habits close most of these gaps. Fit the negative binomial alongside the Poisson every single time, because the comparison costs one line and settles the overdispersion question for you. And keep your syntax in a file rather than clicking through menus, so the next analyst can reproduce the exact analysis.
Frequently Asked Questions
What data can I use with Poisson regression?
Poisson regression needs a non-negative integer count as the outcome: complaints, visits, accidents, defects, absences, or responses such as times eating fast food per week. Every unit must be one whole event, and decimals mean you have a rate or proportion instead. Counts that are naturally bounded, like a 0 to 5 satisfaction scale, are usually better served by ordinal methods.
How do I test for overdispersion in a Poisson model?
Compare the variance of the outcome with its mean, then compute the dispersion statistic from the fitted model as the sum of squared Pearson residuals divided by the residual degrees of freedom. In R that is sum(residuals(m, type = pearson)^2) divided by df.residual(m). In Stata use estat gof. A value near 1 supports the Poisson model; well above 1 means you need quasi-Poisson or negative binomial.
When should I use negative binomial instead of Poisson regression?
Use the negative binomial when dispersion runs clearly above 1, or when the variance grows faster than the mean, which is common with counts above about ten. Its variance function is the mean plus the mean squared times a dispersion parameter kappa, so heavier right tails are handled properly. The likelihood ratio test between the two fits, anova(m_pois, m_nb, test = Chisq), settles close calls.
How do I handle zero-inflated count data?
First compare observed zeros with the zeros your fitted model predicts. A gap of several percentage points suggests a mixture of two processes, and then a zero-inflated Poisson or negative binomial is worth fitting, along with a hurdle model. If the gap is small, keep the standard model, because an unnecessary zero-inflated fit adds parameters and often weakens inference on the predictors.
Which software should I use, and how do I read the output?
R uses glm() with family = poisson, Python uses statsmodels GLM with family = sm.families.Poisson(), SPSS uses Analyze, Generalized Linear Models, Poisson, and Stata uses the poisson command. All four report the coefficient on the log scale. Exponentiate the Estimate column to get incidence rate ratios, and read the Scale or dispersion line to judge fit.
Conclusion
To run Poisson regression for count data properly, do six things in order: confirm the outcome is a non-negative integer count, check missing values and duplicates and fix the reference levels, fit the model with a log link and an offset if exposure varies, compute the dispersion and the residual diagnostics, switch to quasi-Poisson, negative binomial or a zero-inflated model when the evidence points there, and report incidence rate ratios with confidence intervals alongside the fit evidence. Fit the negative binomial alongside the Poisson from the start. It costs one line, and it answers the overdispersion question before a referee asks it.
