Poisson Mixed-Effects Models for Counts

Why Model Repeated Counts With a Poisson Mixed-Effects Model?

Counts arise when each observation records how many times an event occurred in a defined observation window: for example, the number of symptoms reported during a visit or the number of events recorded in a study period. Counts are nonnegative integers, and repeated observations from the same subject may be related. A Poisson mixed-effects model represents that relationship by allowing a subject-specific random effect to influence the expected count.

This tutorial considers counts recorded over comparable, fixed observation windows. The model’s mean is the expected count in that window. A count per unit of exposure requires additional handling; that is a separate modeling issue. The essential distinction here is that the model describes the count distribution conditional on the subject’s random effect. That conditional distribution is the basis for interpreting coefficients and calculating probabilities.

Definition: Let $Y_{ij}$ be the nonnegative integer count for subject $i$ at observation $j$. A Poisson-log mixed-effects model specifies $$ Y_{ij}\mid b_i \sim \operatorname{Poisson}(\lambda_{ij}), \qquad \log(\lambda_{ij}) = \mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta} + \mathbf{z}_{ij}^{\mathsf T}\mathbf{b}_i. $$ Here $\lambda_{ij}>0$ is the subject- and observation-specific conditional mean count, $\mathbf{x}_{ij}$ and $\mathbf{z}_{ij}$ are covariate vectors, $\boldsymbol{\beta}$ is a vector of fixed effects, and $\mathbf{b}_i$ is a vector of subject-level random effects. A standard specification assumes $\mathbf{b}_i\sim N(\mathbf{0},\mathbf{G})$, with independent random effects across subjects. Given the random effects, observations are assumed conditionally independent and to follow the stated Poisson distributions.

The log link ensures that $\lambda_{ij}$ is positive while allowing the linear predictor to include covariates and random effects. For a Poisson variable, its conditional mean and conditional variance are both $\lambda_{ij}$:

E(Y_{ij}\mid b_i)=\lambda_{ij}, \qquad \operatorname{Var}(Y_{ij}\mid b_i)=\lambda_{ij}.

The variance equality is a model assumption, not a universal property of count data. The random effect accounts for a subject-specific shift in the conditional mean; it does not remove the need to assess whether the conditional Poisson distribution is a reasonable description.

Count Probabilities and Conditional Rate Ratios

Once the individual-specific mean $\lambda_{ij}$ is specified, the conditional probability of observing a count $y$ follows from the Poisson probability mass function:

P(Y_{ij}=y\mid b_i) = \frac{\exp(-\lambda_{ij})\lambda_{ij}^{y}}{y!}, \qquad y=0,1,2,\ldots.

Thus, probabilities for a particular subject’s count are obtained by first calculating that subject’s $\lambda_{ij}$, including the subject’s random effect, and then substituting it into the Poisson formula. The probability of zero is $\exp(-\lambda_{ij})$; it is not the same quantity as the expected count.

A fixed-effect coefficient has a multiplicative interpretation on the conditional mean scale. If a covariate $x_k$ increases by one unit while the other covariates and the random effect are held fixed, its conditional mean is multiplied by $\exp(\beta_k)$. For a change of $\Delta$ units, the multiplier is $\exp(\beta_k\Delta)$. This is a conditional rate ratio: it compares expected counts for otherwise identical observations at the same random-effect value. It does not directly describe the change in the probability of any particular count.

Worked Example: Calculating an Individual’s Count Probabilities

Suppose a hypothetical model has an intercept $\beta_0=\log(2)$, a treatment coefficient $\beta_1=\log(1.5)$, and a subject-specific random intercept $b_i=\log(0.8)$. For a treated observation, the linear predictor is

\log(\lambda_{ij})=\log(2)+\log(1.5)+\log(0.8), \qquad \lambda_{ij}=2(1.5)(0.8)=2.4.

For this subject and observation window, the expected count is 2.4. The probability of zero events is

P(Y_{ij}=0\mid b_i)=\exp(-2.4)\frac{2.4^0}{0!} =0.09072.

The probability of exactly one event is

P(Y_{ij}=1\mid b_i)=\exp(-2.4)\frac{2.4^1}{1!} =0.2177.

For exactly three events,

P(Y_{ij}=3\mid b_i) =\exp(-2.4)\frac{2.4^3}{3!} =0.09072\frac{13.824}{6} =0.2090.

These are probabilities conditional on the stated random effect, not probabilities averaged over the subject population. The calculations use a Poisson count over the specified window and the model’s individual-specific mean.

Worked Example: Interpreting a Fixed Treatment Effect

In a hypothetical study, the fitted treatment coefficient is $\hat{\beta}_1=0.4055$, with treatment coded 1 and control coded 0. The conditional treatment-to-control mean ratio is

\exp(\hat{\beta}_1)=\exp(0.4055)\approx 1.500.

At the same random-effect value and the same values of the other covariates, the model therefore estimates a mean count 1.5 times as large under treatment as under control during the observation window. Equivalently, that is a 50.0% higher conditional mean count. It is not a 50-percentage-point increase in the probability of an event, and the model coefficient alone does not establish a causal treatment effect.

For example, if the estimated control mean for a particular random-effect value is 2.0 counts per window, the corresponding treated mean is $2.0(1.500)=3.0$ counts per window. The ratio is a comparison of means, not of the probabilities of observing exactly one, two, or any other particular count.

Likelihood Contributions and Population Averaging

The conditional likelihood contribution for one observation is its Poisson probability mass. For a subject with several observations, conditional independence means that the subject’s conditional likelihood is the product of those probabilities. Taking logarithms gives a convenient sum of log-probabilities:

L_i(\mathbf{b}_i) = \prod_j \frac{\exp(-\lambda_{ij})\lambda_{ij}^{y_{ij}}}{y_{ij}!}, \qquad \ell_i(\mathbf{b}_i) = \sum_j \left[ -\lambda_{ij}+y_{ij}\log(\lambda_{ij})-\log(y_{ij}!) \right].

This is the likelihood conditional on $\mathbf{b}_i$. In a mixed-effects analysis, the random effect is not ordinarily known. The subject’s marginal likelihood contribution is obtained by averaging the conditional likelihood over the assumed random-effects distribution:

L_i = \int L_i(\mathbf{b}_i)\,f(\mathbf{b}_i;\mathbf{G})\,d\mathbf{b}_i,

where $f(\mathbf{b}_i;\mathbf{G})$ is the density of the random effects. The integral explains why treating the random effect as known and integrating over its distribution are different operations. Likewise, a population-averaged count probability is obtained by integrating the individual conditional probability over the random-effects distribution, rather than simply setting the random effect to zero.

Worked Example: A Repeated-Count Likelihood Contribution

Suppose a hypothetical subject has one control observation with count $y_{i1}=1$ and one treated observation with count $y_{i2}=3$. Use $\beta_0=\log(2)$, $\beta_1=\log(1.5)$, and $b_i=\log(0.8)$, as in the previous count-probability calculation. The conditional means are

\lambda_{i1}=2(0.8)=1.6, \qquad \lambda_{i2}=2(1.5)(0.8)=2.4.

The probability of the control count, conditional on this random effect, is

P(Y_{i1}=1\mid b_i) =\exp(-1.6)\frac{1.6}{1!} \approx 0.3230.

The probability of the treated count is

P(Y_{i2}=3\mid b_i) =\exp(-2.4)\frac{2.4^3}{3!} \approx 0.2090.

By conditional independence, the subject’s conditional likelihood contribution is their product:

L_i(b_i)=0.3230(0.2090)\approx 0.06751.

This product is not yet the marginal likelihood contribution for the subject, because it conditions on a particular $b_i$. The marginal contribution would average the product over the random-intercept distribution.

Justification of the Log-Link Interpretation

The rate-ratio interpretation follows directly from the model equation. Consider two observations that differ only in covariate $x_k$ by $\Delta$, with the same random effect and all other covariates held fixed. Their linear predictors differ by $\beta_k\Delta$. Exponentiating each predictor and taking the ratio gives

\frac{\lambda_{\text{new}}}{\lambda_{\text{reference}}} = \frac{\exp(\eta_{\text{reference}}+\beta_k\Delta)} {\exp(\eta_{\text{reference}})} = \exp(\beta_k\Delta),

where $\eta_{\text{reference}}$ is the reference linear predictor. The subject-specific contribution cancels because it is held fixed in this comparison. This derivation also shows why $\beta_k$ itself is a log rate ratio, while $\exp(\beta_k)$ is the rate ratio on the original scale.

Assumptions, Diagnostics, and Limitations

Before using the model’s probabilities or coefficient interpretations, check that its assumptions fit the data structure and scientific question:

Each recorded response is a nonnegative integer count for a clearly defined observation window.
Conditional on the random effects and modeled covariates, each response follows a Poisson distribution with the specified log-linear mean.
Repeated responses are conditionally independent given the subject-level random effects.
The random-effects distribution is appropriately specified, with independent subject-level random effects as assumed by the model.
The covariates and fixed-effects structure adequately represent the systematic changes in the conditional mean.

A useful diagnostic is to compare observed counts with counts simulated from the fitted model, including the repeated observations within subjects. Also examine residual patterns against fitted means, covariates, and observation order. For example, a Pearson residual is $(y_{ij}-\hat{\lambda}_{ij})/\sqrt{\hat{\lambda}_{ij}}$, using a fitted conditional mean. Persistent patterns or a conditional spread larger than the Poisson mean-variance assumption predicts may signal lack of fit. Diagnostics should account for the fact that fitted random effects and other estimated parameters are uncertain; residuals computed using fitted values are not independent observations of model error.

The Poisson model is not suitable merely because the outcome consists of counts. A fixed number of trials with success or failure outcomes may instead call for a binomial model, as discussed in the earlier material on Bernoulli trials versus binomial counts. Nor should a Poisson model be assumed adequate if the conditional variance systematically exceeds the conditional mean, or if the model fails to reproduce important features of the observed counts. Those patterns call for further model assessment rather than an automatic interpretation of the fitted rate ratios.

Common mistake: “The treatment coefficient is 0.4055, so treatment increases the probability of a count by 0.4055.” This confuses a log-mean coefficient with a probability difference. With a log link, the conditional mean ratio is $\exp(0.4055)\approx1.500$. To obtain the probability of a specific count, first calculate the individual-specific mean and then use the Poisson probability mass function.

The model’s main strength is a direct, interpretable link between covariates and the conditional expected count, while subject-level random effects represent repeated-measure dependence. Its limitations follow from its assumptions: the Poisson conditional variance equals its conditional mean, and the chosen random-effects and mean structures must be adequate. A fixed-effect rate ratio is conditional on the random effect and covariates; it is not itself a full description of population-level count probabilities.

Check Your Understanding

Try each question, then compare your reasoning with the solutions below.

If an individual’s conditional mean is 1.2 counts per window, what is the conditional probability of zero counts?
A fixed-effect coefficient is $-0.2231$ for a one-unit covariate increase. What is the corresponding conditional mean ratio, approximately?
Why does multiplying two Poisson probabilities give a repeated subject’s conditional likelihood contribution only when the required conditional independence assumption holds?
Does a conditional treatment mean ratio of 1.5 imply that the probability of exactly one count is 1.5 times as large? Explain.

Solutions. (1) $P(Y=0\mid b)=\exp(-1.2)\approx0.3010$. (2) $\exp(-0.2231)\approx0.8000$, so the conditional mean is about 20% lower per one-unit increase, holding other terms fixed. (3) The product is the joint probability only when the repeated outcomes are conditionally independent; otherwise their joint distribution must account for their dependence. (4) No. The ratio compares conditional expected counts. Probabilities of exactly one count must be calculated from each relevant mean using the Poisson mass function and need not have that same ratio.