Comparing rates and Poisson regression

In the last lesson you turned events and person-time into a single rate with a confidence interval. The interesting epidemiology starts when you have two rates and want to know whether one group is worse off. A TB cohort in Selangor might record one notification rate among people living with HIV and another among HIV-negative members. Do the rates differ, and by how much? This lesson answers that with the rate ratio, then shows how Poisson regression gives the same answer and keeps going when you need to adjust for a confounder.

Recall the notation from 11.1. A rate is events divided by person-time, , often reported per 1000 person-years. We keep for the exposed group and for the unexposed (baseline) group.

The rate ratio as the measure of effect

To compare two rates you divide them. The rate ratio is the rate in the exposed group over the rate in the unexposed group.

A rate ratio of 1 means the two rates are equal. Above 1, the exposed group has the higher rate; below 1, the exposure looks protective. We use the ratio rather than the difference because it is usually more stable across populations and it is what Poisson regression returns. For a rare outcome it also sits close to the risk ratio and the odds ratio.

Confidence interval on the log scale

A rate ratio cannot go below 0 and its sampling distribution is skewed, so you do not build the interval directly on the ratio. You build it on the log scale, where the distribution is close to symmetric, then exponentiate back. The standard error of the log rate ratio depends only on the two event counts, not on the person-time.

From there the 95% interval uses an error factor, the same multiplicative form you would use for a risk ratio or odds ratio.

Only the event counts drive the width. Long follow-up with few events still gives a wide interval, because precision comes from how many events you saw, not how many person-years you banked.

A hypothesis test that two rates are equal

The null hypothesis is that the two rates are the same, which on the log scale means the log rate ratio is 0. A Wald z-test divides the log rate ratio by its standard error.

Compare against the standard normal to read off a P-value. When event counts are small the Wald test is approximate, and you can fall back on an exact comparison of two Poisson counts, which poisson.test provides in R.

Worked example 1: TB rates in two groups

A cohort in the Klang Valley follows two groups for TB notification. Among people living with HIV there are 41 cases in 1200 person-years; among HIV-negative members there are 27 cases in 3100 person-years.

1. Rates: and per 1000 person-years.
2. Rate ratio: .
3. Standard error of the log: .
4. Error factor: , so the interval is to .
5. Test: , well past 1.96, so .

The TB rate among people living with HIV is about 3.9 times the rate in HIV-negative members, with 95% confidence that the true ratio lies between 2.4 and 6.4. That is strong evidence the rates differ.

d1 <- 41; T1 <- 1200
d0 <- 27; T0 <- 3100

rr <- # fill in the rate ratio

rr <- (d1 / T1) / (d0 / T0)
se <- sqrt(1/d1 + 1/d0)
ef <- exp(1.96 * se)
round(c(RR = rr, lower = rr / ef, upper = rr * ef), 3)
#    RR  lower  upper
# 3.923  2.413  6.376
z  <- log(rr) / se
pt <- poisson.test(c(d1, d0), c(T1, T0))   # exact comparison of two rates
round(c(z = z, exact_p = pt$p.value), 4)
# z = 5.515, exact_p < 0.001

import numpy as np
from scipy.stats import norm
d1, T1 = 41, 1200
d0, T0 = 27, 3100

rr = # fill in the rate ratio

rr = (d1 / T1) / (d0 / T0)
se = np.sqrt(1/d1 + 1/d0)
ef = np.exp(1.96 * se)
z  = np.log(rr) / se
p  = 2 * norm.sf(abs(z))
print(round(rr, 3), round(rr / ef, 3), round(rr * ef, 3))
print(round(z, 3), p)
# 3.923 2.413 6.376
# 5.515 ~3.5e-08

The Poisson regression model

The hand calculation works for two groups, but it does not extend cleanly to several exposures or a confounder. Poisson regression does. It is for rates exactly what logistic regression is for odds. The model writes the log rate as a linear predictor.

With one binary exposure , coded 1 for exposed and 0 for baseline, the intercept is the log rate in the baseline group, and is the log rate ratio. The model is additive on the log scale and multiplicative on the rate scale: , baseline rate times rate ratio.

Interpreting a coefficient as a log rate ratio

A Poisson coefficient is a log rate ratio. To read it as a rate ratio, exponentiate it.

The same applies to the confidence interval: fit on the log scale, then exponentiate the two endpoints. For a single binary exposure this returns the identical rate ratio, interval, and P-value you computed by hand, because the model-based standard error of the coefficient is exactly .

The offset term

A regression predicts the count of events, but the count depends on how long each group was followed. We need to model events per unit time. Use and move across:

The term sits in the model with its coefficient fixed at 1. It is not estimated, so it is not a predictor. A term like this is called an offset. It is what converts a model for counts into a model for rates. Forget the offset and you are modelling raw event counts, so a group with more person-time looks higher-risk purely because it was watched longer.

Worked example 2: Poisson regression with an offset

Take the same TB data as two grouped rows, fit the model with log(person-time) as the offset, and read the rate ratio on the exponentiated scale.

events <- c(41, 27)
pyears <- c(1200, 3100)
hiv    <- c(1, 0)          # 1 = HIV positive (exposed), 0 = baseline

m <- glm(events ~ hiv, family = poisson, offset = log(pyears))

m <- glm(events ~ hiv, family = poisson, offset = log(pyears))
exp(cbind(estimate = coef(m), confint.default(m)))
# (Intercept)  baseline rate per person-year (~0.00871)
# hiv          rate ratio 3.923, 95% CI 2.413 to 6.376

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
df = pd.DataFrame({"events": [41, 27],
                   "pyears": [1200, 3100],
                   "hiv":    [1, 0]})

m = smf.glm("events ~ hiv", data=df, family=sm.families.Poisson(),
            offset=np.log(df["pyears"])).fit()

m = smf.glm("events ~ hiv", data=df, family=sm.families.Poisson(),
            offset=np.log(df["pyears"])).fit()
print(np.exp(m.params))
print(np.exp(m.conf_int()))
# hiv: rate ratio 3.923, 95% CI 2.413 to 6.376

Adjusting for a confounder

A confounder is linked to both the exposure and the outcome and is not on the causal path between them. In the TB cohort, age could be one: older members may have both higher TB rates and a different HIV mix. You adjust by adding the confounder as another term in the model. The coefficient on the exposure then becomes the rate ratio holding the confounder fixed.

events <- c(30, 11, 20, 7)
pyears <- c(700, 500, 1500, 1600)
hiv    <- c(1, 1, 0, 0)
older  <- c(1, 0, 1, 0)   # 1 = aged 35+, 0 = under 35

m2 <- glm(events ~ hiv + older, family = poisson, offset = log(pyears))

m2 <- glm(events ~ hiv + older, family = poisson, offset = log(pyears))
round(exp(cbind(RR = coef(m2), confint.default(m2))), 3)
# hiv row   = rate ratio for HIV, adjusted for age group
# older row = rate ratio for age, adjusted for HIV

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
df2 = pd.DataFrame({"events": [30, 11, 20, 7],
                    "pyears": [700, 500, 1500, 1600],
                    "hiv":    [1, 1, 0, 0],
                    "older":  [1, 0, 1, 0]})

m2 = smf.glm("events ~ hiv + older", data=df2, family=sm.families.Poisson(),
             offset=np.log(df2["pyears"])).fit()

m2 = smf.glm("events ~ hiv + older", data=df2, family=sm.families.Poisson(),
             offset=np.log(df2["pyears"])).fit()
print(np.exp(m2.params))
print(np.exp(m2.conf_int()))
# hiv = rate ratio for HIV adjusted for age group

Read the rate ratio for HIV from the hiv row and compare it to the crude 3.92. If it barely moves, age was not carrying much of the association and is not an important confounder here. If it shifts a lot, the crude figure was mixing the HIV effect with the age effect, and the adjusted value is the one to report. The model assumes the HIV rate ratio is the same in both age groups; checking that assumption (effect modification) is a later topic.

When Poisson is inadequate: overdispersion

The Poisson model assumes the variance of the count equals its mean. Real data often have more spread than that, because events cluster (households, outbreaks, repeat episodes) or because an important covariate is missing. This is overdispersion. The symptom is a residual deviance far larger than the residual degrees of freedom, or a dispersion statistic above 1. The danger is that the rate ratios stay roughly right but the standard errors come out too small, so confidence intervals are too narrow and P-values too optimistic.

Two standard fixes keep the same rate-ratio interpretation. Quasi-Poisson multiplies the standard errors by the square root of the estimated dispersion, widening the intervals. Negative binomial regression adds a separate variance parameter and models the extra spread directly. Both are one line away from the Poisson fit.

events <- c(30, 11, 20, 7)
pyears <- c(700, 500, 1500, 1600)
hiv    <- c(1, 1, 0, 0)
older  <- c(1, 0, 1, 0)

m_qp <- glm(events ~ hiv + older, family = quasipoisson, offset = log(pyears))

m_qp <- glm(events ~ hiv + older, family = quasipoisson, offset = log(pyears))
round(summary(m_qp)$dispersion, 3)   # > 1 would signal overdispersion

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
df2 = pd.DataFrame({"events": [30, 11, 20, 7],
                    "pyears": [700, 500, 1500, 1600],
                    "hiv":    [1, 1, 0, 0],
                    "older":  [1, 0, 1, 0]})

m_nb = smf.glm("events ~ hiv + older",
               family=sm.families.NegativeBinomial(alpha=1.0), ...).fit()

m_nb = smf.glm("events ~ hiv + older", data=df2,
               family=sm.families.NegativeBinomial(alpha=1.0),
               offset=np.log(df2["pyears"])).fit()
print(round(np.exp(m_nb.params["hiv"]), 3))   # rate ratio under negative binomial

Common mistakes

• Building the rate-ratio interval on the ratio scale. The sampling distribution of a rate ratio is skewed. Add and subtract on the ratio scale and the lower limit can fall below 0. Always work on the log scale and exponentiate.
• Dropping the offset, or passing raw person-time. Without the model fits counts, not rates, and a group followed longer looks higher-risk for the wrong reason. The offset is the log of person-time, not person-time itself.
• Reading the raw coefficient as a rate ratio. A Poisson coefficient is a log rate ratio. A value of 0.53 is not a rate ratio of 0.53; it is . Exponentiate before you interpret.
• Trusting narrow intervals under overdispersion. If the dispersion statistic is well above 1, the Poisson standard errors are too small. Refit with quasi-Poisson or negative binomial before quoting a confidence interval.
• Confusing precision from person-time with precision from events. The standard error of the log rate ratio depends only on the event counts. More person-years with no extra events does not tighten the interval.

Tips

• Report rate ratios with their confidence intervals, not P-values alone. The interval shows both the size of the effect and how well the data pin it down.
• Code the exposure so the baseline group is 0. The intercept is then the log baseline rate, and each other coefficient is a clean log rate ratio against it.
• Fit the crude model first, then add the confounder. The shift in the rate ratio tells you how much confounding the variable was carrying.
• For a quick check on published tables, remember the limits are multiplicative: rate ratio divided by lower limit should equal upper limit divided by rate ratio.
• When counts are clustered (households, clinics, repeat episodes), expect overdispersion and reach for negative binomial or a quasi-Poisson scale rather than plain Poisson.