How do you model counts and rates?

Start here

What makes a count different from a measurement?

A count is the number of times something happened: seizures in a month, infections in an ICU, admissions in a year. It can only be 0, 1, 2 and so on. It cannot be negative, and when the average is small most values pile up near 0.

Three features of counts cause problems for the t-test and linear regression. The distribution is skewed when the mean is low. The variance rises with the mean, so groups with more events vary more. And the count depends on how long, or how many people, you watched: a patient followed for 6 months has more chance to record seizures than a patient followed for 3.

The Poisson distribution

When does a count follow a Poisson distribution?

When there are many opportunities for an event, each with a small probability, the events happen independently of each other, and the underlying rate stays the same, the number of events follows a Poisson distribution.

P(k events) = e^−μ × μ^k ÷ k!

k

number of events: 0, 1, 2, …

μ

mean number of events in the interval (in the figure, 2 deaths a year)

k!

k factorial: 1 × 2 × … × k, with 0! = 1

mean = variance

both equal μ, so the SD is √μ

Worked example

With a mean of 2 events a year:

• P(0) = e⁻² = 0.1353. In about 14 years out of 100 there would be no such death.
• P(1) = e⁻² × 2 ÷ 1 = 0.2707. P(2) = e⁻² × 4 ÷ 2 = 0.2707.
• P(5 or more) = 1 − [P(0) + P(1) + P(2) + P(3) + P(4)] = 1 − 0.9473 = 0.0527.

Uses of the Poisson distribution

• Rare events in large populations: cases of a disease per year in a district, deaths per 1,000 person-years.
• Counting in the laboratory and in nuclear medicine: counts from radioactive decay, colonies on a plate, cells in a counting chamber. The SD of a count N is √N. 10,000 counts have an SD of 100 (1%), and 100 counts an SD of 10 (10%).
• Checking whether events cluster: if events occur more in bunches than a Poisson distribution predicts, the variance will be larger than the mean.
• Poisson regression: modelling counts and rates with several predictors, described below.

For a large mean the Poisson distribution becomes close to a normal distribution with mean μ and SD √μ. With μ = 16, the exact P(10 or fewer) is 0.0774 and the normal approximation with continuity correction gives 0.0846.

Interactive

Change the average seizure count and watch the Poisson probabilities and the variance change

In a Poisson distribution the variance equals the mean, so the spread rises as the mean rises.

Rates and the offset

Why compare events per person-time and not the counts?

A count depends on how long, and on how many people, you observed. A rate divides the count by the time at risk: infections per 1,000 catheter-days, deaths per 1,000 person-years, seizures per patient-month.

Poisson regression with an offset

Poisson regression models the logarithm of the expected count. To model a rate, the logarithm of each unit's time at risk is added to the model with its coefficient fixed at 1. This term is the offset. The SPSS manual describes it as a “structural” predictor whose coefficient is not estimated but assumed to be 1.

ln(expected count) = ln(time at risk) + b₀ + b₁x₁ + … + bₖxₖ

which is the same as   ln(expected count ÷ time at risk) = b₀ + b₁x₁ + … + bₖxₖ

ln(time at risk)

the offset: catheter-days, person-years, months of follow-up, population × years

e^b₀

the rate when every predictor is 0

e^bᵢ

the incidence rate ratio (IRR) for a 1-unit rise in xᵢ, adjusted for the other predictors

Worked example: two ICUs

IRR = (18 ÷ 3,600) ÷ (9 ÷ 1,200) = 5.0 ÷ 7.5 = 0.67

SE of ln IRR = √(1/18 + 1/9) = 0.408   95% CI = e^{ln 0.667 ± 1.96 × 0.408} = 0.30 to 1.48

p = 0.32. A Poisson regression with ICU as the only predictor and ln(catheter-days) as the offset gives the same IRR, confidence interval and p-value. Without the offset, the ratio of counts (2.0) would point the wrong way.

Incidence, prevalence and person-time are introduced in the risk and survival chapter.

Interactive

Change infections and catheter-days and watch the counts and the rates disagree

The incidence rate ratio uses the time at risk. Its confidence interval depends on the number of events.

Reading the output

How do you read an incidence rate ratio?

An incidence rate ratio compares event rates. An IRR of 0.78 means the rate in one group is 78% of the rate in the reference group, a 22% lower rate, after adjusting for the other predictors.

• Per category: each category is compared with the reference category, as in logistic regression.
• Per unit: for a continuous predictor the IRR is per 1-unit rise. For k units use e^{k × b}.
• Rate, not risk: the IRR compares events per unit of time. It allows for repeated events in the same patient and for different follow-up.
• Absolute numbers: a ratio does not show how many events are prevented. Give the rates in each group as well, for example seizures per patient-month (STROBE item 16c asks authors to consider translating relative risks into absolute risks).

Poisson regression can also estimate risk ratios for a yes/no outcome. Zou (2004) proposed Poisson regression with robust standard errors for prospective studies with binary data, the “modified Poisson” approach, because the usual Poisson standard errors are wrong for a 0/1 outcome. The logistic regression chapter explains when a risk ratio is needed instead of an odds ratio.

Overdispersion

What if the counts vary more than a Poisson distribution allows?

Poisson regression assumes that patients with the same predictor values have the same underlying rate. Real patients differ: some have frequent seizures and some rare ones. Mixing patients with different rates makes the variance larger than the mean. This is overdispersion.

Why it matters

Poisson regression calculates its standard errors from the assumption that variance equals mean. When the variance is larger, the standard errors are too small, the confidence intervals too narrow and the p-values too small. The rate ratio itself is usually estimated reasonably, but its precision is overstated.

How to detect it

1. Variance ÷ mean of the counts in each group. In the teaching trial below: 34.4 ÷ 7.6 = 4.5 with standard treatment and 72.6 ÷ 6.17 = 11.8 with the add-on drug. (With unequal follow-up this is only a rough check.)
2. Pearson χ² ÷ residual degrees of freedom after fitting Poisson regression. Close to 1 is consistent with a Poisson distribution. The teaching trial gives 813.07 ÷ 118 = 6.89. The deviance ÷ df, 5.48, is a similar check.
3. Likelihood ratio test of α = 0: compares the negative binomial model with the Poisson model. Teaching trial: χ² = 332.4, p < 0.001.
4. Observed and expected frequencies: compare how many patients had 0, 1, 2 … events with the numbers each model predicts.

The fix: negative binomial regression

Negative binomial regression uses the same log link and offset, and gives incidence rate ratios in the same way. It adds one parameter, α, for the extra variation between patients:

variance = μ + α × μ²

When α = 0 the model is Poisson regression. Two other corrections keep the Poisson estimates but enlarge the standard errors: quasi-Poisson regression, which multiplies them by √(Pearson χ² ÷ df), and robust standard errors. In the teaching trial all three corrections give similar intervals, all wider than the Poisson interval.

Interactive

Increase the extra variation and watch the Poisson confidence intervals stop covering the true rate ratio

A 95% confidence interval should contain the true value in about 95 of 100 studies. With overdispersion, Poisson intervals do so less often.

Zeros

What if there are more zeros than the model expects?

Sometimes the zeros come from two sources. Some people cannot have the event at all, such as non-smokers when the outcome is cigarettes smoked. Others could have it but did not, by chance. A zero-inflated model treats these two sources separately.

• Zero-inflated model: Lambert (1992) described zero-inflated Poisson regression: with probability p the only possible observation is 0, and with probability 1 − p a Poisson count is observed. Both parts can depend on predictors. A zero-inflated negative binomial model uses a negative binomial count part.
• Hurdle model: one part models whether the count is zero or not, and a second part models the size of the count among those above zero.
• Check the negative binomial model first. Many zeros are not the same as excess zeros. In the seizure trial the negative binomial model predicts 10.0 patients with no seizures against 9 observed, so no zero-inflated model is needed.
• Choose a zero-inflated model only when there is a clinical reason for a group that cannot have the event.

Multi-way tables

How do you analyse a table with three categorical variables?

A log-linear model is Poisson regression applied to the cell counts of a contingency table. There is no single outcome variable. The interaction terms describe which variables are associated.

The models, from simplest to saturated

Write hospital as H, technique as T and infection as O. The notation lists the highest interactions included: [HT][HO] means the model contains H, T, O, the H × T term and the H × O term.

G² (the deviance) compares each model's expected counts with the observed counts. A small G² with a large p-value means the model fits. The conditional independence model [HT][HO] fits, and adding the T × O term improves G² by only 0.05 (p = 0.82). The technique–infection association in the pooled table (G² = 6.89, 1 df, p = 0.009) disappears once hospital is in the model.

ln(expected count) = λ + λ^H + λ^T + λ^O + λ^HT + λ^HO

Decide

Which count model fits your data?

Answer the questions in order. Each answer links to its section.

The formula

Worked example: the seizure trial

Made-up teaching data: the trial at the top of this chapter. 60 patients per group, followed for 3 to 6 months. Fitted in Python with statsmodels, with ln(months of follow-up) as the offset.

1. Crude rate ratio: (370 ÷ 282) ÷ (456 ÷ 271) = 1.312 ÷ 1.683 = 0.78.
2. Poisson regression: B = −0.2488, SE 0.070. IRR = e^−0.2488 = 0.78, 95% CI e^{−0.2488 ± 1.96 × 0.070} = 0.68 to 0.89, p < 0.001.
3. Check dispersion: Pearson χ² 813.07 on 118 df, ratio 6.89. The counts vary about seven times more than a Poisson distribution allows.
4. Negative binomial regression: B = −0.2422, SE 0.165, α = 0.67. IRR = 0.78, 95% CI 0.57 to 1.09, p = 0.14.
5. Compare the models: likelihood ratio test of α = 0, χ² = 332.4, p < 0.001. AIC 1041.3 (Poisson) against 710.9 (negative binomial).
6. Other corrections: quasi-Poisson 95% CI 0.54 to 1.12 (p = 0.18). Poisson with robust standard errors 0.53 to 1.15 (p = 0.21).
7. Zeros: 9 observed, 0.4 expected under Poisson, 10.0 under the negative binomial model. No zero-inflated model is needed.
8. Absolute difference: 1.683 − 1.312 = 0.37 fewer seizures per patient-month, 37.1 per 100 patient-months.

Conclusion: the point estimate is the same, a 22% lower seizure rate, but once the variation between patients is allowed for, the data are compatible with anything from a 43% reduction to a 9% increase.

In a real paper

Count models in two published cohort studies

How to run it

Running Poisson, negative binomial and log-linear models in SPSS, R and jamovi

Menu paths are from the IBM SPSS Statistics 27 manuals. Functions and arguments are from the CRAN reference manuals and the R documentation.

SPSS

1. Create the offset: Transform > Compute Variable, with ln_months = LN(months). The offset must be the logarithm of the time at risk.
2. Analyze > Generalized Linear Models > Generalized Linear Models. On the Type of Model tab, under Counts, choose Poisson loglinear.
3. On the Response tab choose the count. On the Predictors tab add factors and covariates and put ln_months in Offset. Add the effects on the Model tab.
4. The goodness-of-fit statistics include the deviance and Pearson chi-square. Divide Pearson chi-square by the residual degrees of freedom to check dispersion. The IRR is e^B.
5. For negative binomial regression, the Negative binomial with log link option fixes the ancillary parameter at 1. To estimate it, choose Custom, distribution Negative binomial, link Log, and Estimate value for the parameter.
6. Log-linear models: Analyze > Loglinear > Model Selection fits hierarchical log-linear models with backward elimination. Analyze > Loglinear > General fits a specified model.

R

fit_p <- glm(seizures ~ drug + offset(log(months)), family = poisson, data = d)
exp(cbind(IRR = coef(fit_p), confint.default(fit_p)))
sum(residuals(fit_p, type = "pearson")^2) / df.residual(fit_p)   # dispersion ratio
performance::check_overdispersion(fit_p)

fit_nb <- MASS::glm.nb(seizures ~ drug + offset(log(months)), data = d)   # estimates theta = 1/alpha
exp(cbind(IRR = coef(fit_nb), confint.default(fit_nb)))
fit_q  <- glm(seizures ~ drug + offset(log(months)), family = quasipoisson, data = d)
lmtest::coeftest(fit_p, vcov. = sandwich::vcovHC(fit_p, type = "HC0"))    # robust standard errors

pscl::zeroinfl(cigs ~ age | age, dist = "negbin", data = z)   # count part | zero part
performance::check_zeroinflation(fit_nb)

m1 <- glm(n ~ H*T + H*O, family = poisson, data = cells)       # log-linear: [HT][HO]
m2 <- update(m1, . ~ . + T:O)
anova(m1, m2, test = "Chisq")

jamovi

1. Log-linear models: Analyses → Frequencies → Log-Linear Regression, with the categorical variables as factors and the cell counts as the counts variable.
2. Poisson and negative binomial regression need the GAMLj module, installed from the jamovi library. It adds a Linear Models menu containing Generalized Linear Models, with Poisson, Poisson (overdispersion) and Negative binomial model types.
3. GAMLj's Offset variable enters the model with its coefficient fixed at 1, so supply ln(time at risk), not the time itself.

How to report it

What to write in the thesis and the paper

Give the counts, the time at risk and the rates, say how you checked dispersion, then give the rate ratio with its confidence interval.

Methods

• The count outcome and how time at risk was measured (person-years, catheter-days, months of follow-up), and that its logarithm was the offset.
• The model and why: Poisson, negative binomial, quasi-Poisson or robust standard errors, and how overdispersion and excess zeros were assessed (STROBE item 12a: describe all statistical methods).
• Predictors adjusted for and why (STROBE item 16a), reference categories, software and version.

Results

• Events, time at risk and rate in each group, then unadjusted and adjusted IRRs with 95% CIs.
• The dispersion statistic or test result, and absolute rate differences where they help (STROBE item 16c).

A worked sentence using this chapter's teaching numbers

Over 553 patient-months of follow-up there were 456 seizures with standard treatment (1.68 per patient-month) and 370 with the add-on drug (1.31 per patient-month). Seizure counts were overdispersed (Pearson χ² ÷ df 6.89; likelihood ratio test of α = 0, p < 0.001), so negative binomial regression with ln(months of follow-up) as the offset was used. The add-on drug was associated with a 22% lower seizure rate (incidence rate ratio 0.78, 95% CI 0.57 to 1.09, p = 0.14), or 0.37 fewer seizures per patient-month.

By specialty

Count model problems by specialty

Common myths

Common mistaken beliefs

• ✗Myth: More events in one group means a higher rate.

  ✓Only if time at risk is the same. ICU A had twice ICU B's infections and a lower rate: 5.0 against 7.5 per 1,000 catheter-days.

• ✗Myth: A very small p-value from Poisson regression settles the question.

  ✓Check dispersion first. In the teaching trial Poisson regression gave p < 0.001 and negative binomial regression p = 0.14 for the same rate ratio of 0.78.

• ✗Myth: Many zeros always need a zero-inflated model.

  ✓A negative binomial model often predicts them. In the seizure trial it expected 10.0 zeros against 9 observed.

• ✗Myth: The offset is one more predictor to adjust for.

  ✓Its coefficient is fixed at 1, and it must be the logarithm of the time at risk.

• ✗Myth: An incidence rate ratio is the same as a risk ratio.

  ✓The IRR compares events per unit of person-time and counts every event. The risk ratio compares the proportion of people with at least one event.

• ✗Myth: Counts can be analysed with a t-test on the mean count.

  ✓A t-test ignores unequal follow-up and the link between the mean and the variance, and gives no rate ratio.

Why this matters to you

Five consequences of getting count models wrong

1. Protocol sent back

   A plan that compares "number of infections" between wards without catheter-days or patient-days as the denominator.

2. A wrong conclusion in the thesis

   Poisson regression on overdispersed seizure counts gives a p-value that is far too small, and a drug is declared effective on weak evidence.

3. Questions you cannot answer at the viva

   "How did you allow for different follow-up?" and "Did you check whether the variance was larger than the mean?" need the offset and a dispersion statistic.

4. Misreading papers

   Accepting a rate ratio with no statement of the time at risk, or confidence intervals that are narrower than the number of events allows.

5. Harm to patients

   An infection-control unit that treats more catheter-days as "more infections" may change a practice that was working.

History

History of count models

1. 1837

   Siméon-Denis Poisson

   Published a monograph on probabilities in which he derived the distribution now named after him, as a limit of the negative binomial distribution.

2. 1898

   Ladislaus von Bortkiewicz

   Published data on deaths from horse kicks in corps of the Prussian army, each observed over 20 years, in his "law of small numbers". Armitage later noted that the near-Poisson overall distribution could hide considerable variation between corps.

3. 1920

   Major Greenwood and G. Udny Yule

   Published "An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents".

4. 1972

   John Nelder and Robert Wedderburn

   Published "Generalized linear models", which fitted normal, binomial, Poisson (contingency tables) and gamma models by one method.

5. 1989

   C. B. Dean and Jerald Lawless

   Published "Tests for detecting overdispersion in Poisson regression models" in the Journal of the American Statistical Association.

6. 1992

   Diane Lambert

   Described zero-inflated Poisson regression for count data with excess zeros, with an application to soldering defects on printed wiring boards.

7. 2004

   Guangyong Zou

   Proposed modified Poisson regression, with robust standard errors, to estimate relative risks directly in prospective studies with binary data.

Guidelines

Count models in reporting guidelines

Reading a paper

Six checks on any count model

• Is the time at risk (person-years, catheter-days, follow-up) stated, and was its logarithm used as an offset when it differed between units?
• Are the event counts and rates given for each group, not only the rate ratio?
• Was overdispersion checked, and was negative binomial regression, quasi-Poisson or robust standard errors used when it was present?
• If there are many zeros, were they compared with the zeros the model predicts, and is there a clinical reason for any zero-inflated model?
• Are the confidence intervals plausible for the number of events?
• For a multi-way table, was the association checked within strata as well as in the pooled table?

Why is it called that?

Word surgery

Exam answers

Model answer outlines

Past papers

TNMGRMU questions on poisson and log-linear models

Questions from MD, MS, DM, MCh and Ph.D. papers that this chapter helps you answer. Tap a sitting to open the university's own question paper.

1. ×6Poisson distributionMD Nuclear Medicine · May 2023MD Nuclear Medicine · May 2022MD Nuclear Medicine · Jul 2021MD Nuclear Medicine · May 2019MD Community Medicine · Apr 2016MD Community Medicine · Sep 1992