# Pair-matched count regression in R with offset?

I need to conduct a pair-matched count regression in R.

Specifically, my data looks as follows:

ID Exposure Occurrence Area
1 0 7 50
1 1 5 55
2 0 6 70
2 1 8 74
3 0 4 45
3 1 9 62

I need the outcome to be based on a relative rate per square feet (area), so I want to use an offset.

While I would think this wouldn't be too difficult, I am hitting a wall on my Google search as pair-matched data is so frequently done with a binary outcome.

Do you know a good package or a way to be able to do this?

You should be able to do this with a mixed model with a count response (e.g. Poisson or negative binomial): you want the "standard" count-GLM-with-offset model with random variation in the intercept across IDs:

$$\begin{split} \eta_{ij} & = \beta_0 + b_i + \beta_1 E_{ij} + \log(A_{ij}) \\ b_i & \sim N(0, \sigma^2_b) \\ O_{ij} & \sim \textrm{Poisson}(\exp(\eta_{ij})) \end{split}$$ where $$i$$ indexes IDs and $$j = \{1,2\}$$ indexes observations within IDs. $$\beta_0$$ denotes the overall (population-level) intercept; $$b_i$$ denotes the random effects, i.e. the deviation of the intercept from the population-level value for each group. (Most mixed-model software will let you estimate the conditional modes of the $$b_i$$ distributions, equivalent to BLUPs (best linear unbiased predictors) in the linear mixed model case.) $$A_{ij}$$ is the area for the $$i,j$$th observation (i.e. observation $$j$$ within group $$i$$), $$O_{ij}$$ is the occurrence (count response).

If you specified that the model be fitted by restricted maximum likelihood (which is possible for GLMMs in glmmTMB and possibly some other packages), this specification is exactly analogous to a paired t-test, but with Poisson rather than Gaussian responses. If you use maximum likelihood or Bayesian methods (which are more common), it's still pretty close to a "paired Poisson t-test" equivalent.

Since you asked how to do this in R: in lme4 (for example) it would be

glmer(observation ~ exposure + offset(log(area)) + (1|ID),
family = poisson,
data = ...)


Similar models can be fitted lots of different packages/functions, some with very similar interfaces (lme4::glmer, glmmTMB::glmmTMB), some with different interfaces (GLMMadaptive), some Bayesian (MCMCglmm, rstanarm, brms), etc.

• Hi, can you clarify what some of the variables above are standing in for? For example, bi compared to B0? Commented Apr 22, 2022 at 19:06
• @revere2323 $\beta_0$ is the overall model intercept (i.e., the fixed intercept), and $b_i$ is the randomly varying effect for each unit $i$. $b_i$ is specified to have a Normal distribution, and the fact that it varies across individuals and has a probability distribution associated with it is what makes this a mixed model. I infer that $O_{ij}$ is the occurrence for observation $j$ of unit $i$, and $E_{ij}$ is the exposure. $\beta_1$ is the coefficient on exposure, i.e., the ID-conditional exposure effect. This is all described in the paragraph under the regression equation.
– Noah
Commented Apr 22, 2022 at 20:09
• Thanks, Noah. Some of what's in your comment was edited in shortly before you commented . Commented Apr 22, 2022 at 21:22