Fitting Poisson GLMMs with glmer

Why Fit a Poisson Random-Intercept Model?

Event counts collected repeatedly within clinics, schools, or other clusters are usually related within a cluster. A clinic that tends to have higher counts in one period may also tend to have higher counts in other periods. A Poisson generalized linear mixed model (GLMM) represents this clustering with a random intercept, while allowing measured predictors to explain systematic differences in expected counts.

The focus here is fitting and reading a Poisson random-intercept model with lme4::glmer. The example uses event counts per clinic-period, with equal observation periods; it does not model differing exposure times. A random intercept accounts for cluster-level variation, but it does not by itself establish that a predictor caused a change in counts.

Model and Assumptions

Definition: For count $Y_{ij}$ from observation $i$ in cluster $j$, a Poisson random-intercept model specifies $$ Y_{ij}\mid b_j \sim \operatorname{Poisson}(\lambda_{ij}), \qquad \log(\lambda_{ij})=\beta_0+\beta_1 x_{ij}+b_j,\qquad b_j\sim N(0,\tau^2). $$ Here $x_{ij}$ is a predictor, $\beta_0$ and $\beta_1$ are fixed-effect coefficients, and $b_j$ is the cluster-specific random intercept. Conditional on the random effects and predictors, observations are assumed independent. The conditional Poisson mean and variance are both $\lambda_{ij}$.

The model requires nonnegative integer outcomes and an appropriate conditional mean structure. The random intercepts are assumed to follow a normal distribution with mean zero and variance $\tau^2$, independently across clusters, and to be independent of the predictors under the model. There must be multiple clusters to estimate between-cluster variation, and predictors must be coded as intended. In the example, every clinic-period has the same duration, so a count mean is comparable across rows without an exposure adjustment.

These are modeling assumptions, not consequences of fitting the software. In particular, the Poisson variance assumption applies conditional on the random effect. A fitted random intercept can explain dependence and extra variation due to persistent differences between clinics, but it does not guarantee that all remaining conditional variation is Poisson.

Worked Example: Specify and Fit the Model

Suppose a hypothetical health network records the number of acute events in each of four equal-length clinic-periods at 30 clinics. Each clinic contributes periods under routine care and under an enhanced program. Let events be the count, treatment identify the program, and clinic identify the cluster. Before fitting, confirm that counts are integers at least zero, that each row represents one clinic-period, and that the clinic identifiers correctly group repeated observations.

Set routine care as the reference category and fit the following model. The displayed data-frame dimensions describe the hypothetical full dataset; the output numbers below are an illustrative software result, not values calculated from a shortened data excerpt.

\log(\lambda_{ij})=\beta_0+\beta_1 I(\text{enhanced}_{ij})+b_j

treatment is a factor, so explicitly setting its reference level makes the intercept refer to routine care. The model can be fitted as follows:

dat$treatment <- relevel(factor(dat$treatment), ref = "routine")
fit <- lme4::glmer(events ~ treatment + (1 | clinic),
family = poisson(link = "log"), data = dat)
summary(fit)

The term (1 | clinic) requests a clinic-specific random intercept. With glmer's default settings, the model is fit by maximum likelihood using a Laplace approximation to the likelihood integral. A plausible illustrative fixed-effect table is:

Term	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	1.099	0.080	13.74	< 0.001
treatmentenhanced	-0.223	0.100	-2.23	0.026

Suppose the reported standard deviation of the clinic random intercept is 0.400. It is measured on the log-mean scale, so its variance is $0.400^2=0.160$. The $z$ statistic for treatment is the estimate divided by its standard error: $-0.223/0.100=-2.23$. The two-sided Wald $p$-value is approximately 0.026 under a standard normal reference distribution.

Interpret the Fixed Effects on the Count Scale

The intercept is the conditional log expected count for the reference treatment when the random intercept is zero. Exponentiating converts a log mean to a mean count. The treatment coefficient is a conditional log mean ratio: exponentiating it gives the expected-count ratio comparing enhanced with routine care at the same value of the clinic random intercept.

Worked Example: Convert the Treatment Estimate to a Ratio

For the illustrative fitted table, the intercept is 1.099 and the enhanced-treatment coefficient is $-0.223$. At a clinic with $b_j=0$, the estimated routine-care mean is:

\hat{\lambda}_{ij}=\exp(1.099)\approx 3.001

This is approximately 3.001 events per clinic-period. The enhanced-to-routine conditional mean ratio is:

\exp(-0.223)\approx 0.800

Thus, at the same clinic random intercept, the fitted enhanced-program mean is about 0.800 times the routine-care mean, or 20.0% lower. At $b_j=0$, this gives $\exp(1.099-0.223)=\exp(0.876)\approx 2.401$ events per clinic-period. The estimate describes an adjusted association in the fitted model; it is not, by itself, evidence of a causal program effect.

Worked Example: Predict for a Clinic with a Nonzero Random Intercept

Suppose a particular clinic has an estimated random intercept of $b_j=0.200$. This shifts both treatment-specific log means upward by 0.200. The fitted routine-care mean for this clinic is:

\exp(1.099+0.200)=\exp(1.299)\approx 3.666

Equivalently, $\exp(1.099)\exp(0.200)\approx 3.001(1.2214)\approx 3.666$ events per clinic-period. The enhanced mean is $\exp(1.099-0.223+0.200)=\exp(1.076)\approx 2.933$ events per clinic-period. The ratio remains $\exp(-0.223)\approx0.800$, because the same random intercept enters both means and cancels in their ratio. A clinic-specific prediction uses an estimated clinic effect; a prediction for a new clinic cannot use a previously estimated $b_j$.

Uncertainty, Marginal Means, and a Short Justification

A common summary reports a Wald confidence interval for a coefficient as its estimate plus or minus 1.96 standard errors, then exponentiates the endpoints for a ratio interval. For the treatment coefficient, the approximate 95% interval on the log scale is $-0.223\pm1.96(0.100)$, or $[-0.419,-0.027]$. Exponentiating gives a conditional mean-ratio interval of approximately $[0.658,0.973]$. Under the model and the Wald approximation, this interval is compatible with conditional expected counts ranging from about 2.7% to 34.2% lower under enhanced care. The interval procedure has 95% coverage in repeated samples under its assumptions; it does not assign a 95% probability to this particular interval containing the fixed parameter.

The coefficient table's intercept describes the cluster with $b_j=0$, not the average mean across the random-effect distribution. This distinction follows directly by integrating over the normal random intercept:

E(Y_{ij}\mid x_{ij})=E_{b_j}\{\exp(\beta_0+\beta_1x_{ij}+b_j)\} =\exp(\beta_0+\beta_1x_{ij})E\{\exp(b_j)\} =\exp(\beta_0+\beta_1x_{ij}+\tau^2/2).

The last equality uses the moment-generating function of a normal random variable: if $b_j\sim N(0,\tau^2)$, then $E\{\exp(b_j)\}=\exp(\tau^2/2)$. For the illustrative fit, the population-average routine mean is therefore approximately $3.001\exp(0.160/2)=3.001\exp(0.080)\approx3.251$ events per clinic-period. This is larger than the mean at $b_j=0$, because exponentiation is nonlinear. The treatment ratio remains the same here because the same random-intercept distribution is assumed for both treatment conditions.

Diagnostics, Alternatives, and Limitations

A successful optimizer message is necessary but not sufficient evidence that the model is suitable. Check that the fitted model corresponds to the intended observation unit and cluster structure, review convergence warnings, and inspect whether the estimated random-effect variance is near zero. A singular-fit warning can indicate that the data provide little evidence for the specified random-effect variation; it can also occur when the model is too ambitious for the number or information content of clusters.

Compare observed and fitted counts across clinics, and check whether the model systematically misses many zeros, the upper tail, or particular clusters. Simulation-based residual diagnostics can help assess conditional fit. Residual patterns, strong unexplained heterogeneity, or variance beyond what the fitted model captures are reasons to reconsider the specification. Diagnostics do not prove the model assumptions; they can reveal incompatibility. Also check that there are enough independent clusters for reliable estimation and that fitted predictions are not being interpreted beyond the observed design.

A Poisson regression without a random effect is simpler, but treating clustered observations as independent can understate uncertainty when within-cluster dependence is present. A Poisson random-intercept model is useful when persistent cluster-level differences are a plausible source of dependence. If substantial conditional overdispersion remains, a different count distribution or structure may be needed; adding a random intercept is not a universal remedy. The random-intercept model also assumes that clusters differ through a common shift on the log-mean scale. It does not allow predictor effects to vary across clinics, and it can be sensitive to an incorrect random-effect distribution or an inadequate mean structure.

Common pitfall: Reading $\beta_1=-0.223$ as “0.223 fewer events” is incorrect. The coefficient is on the log-mean scale. Exponentiate it to obtain the conditional mean ratio, $\exp(-0.223)\approx0.800$, and describe that ratio in context. Also, do not call the intercept's exponentiated value the overall average count: it is the mean at $b_j=0$, while the population-average mean integrates over the random-intercept distribution.

Check Your Understanding

Use the model and illustrative coefficient table above. Suggested solutions are included after the questions.

What does the treatment coefficient estimate, and why is it not a difference of 0.223 events?
Calculate the approximate 95% Wald interval for the treatment mean ratio using the reported estimate and standard error.
For a clinic with $b_j=0.200$, calculate the routine-care fitted mean to three significant figures.
Why is the marginal routine-care mean not equal to $\exp(1.099)$?

Solutions. (1) It estimates the log of the enhanced-to-routine conditional expected-count ratio; exponentiating gives $\exp(-0.223)\approx0.800$, so it is a ratio, not a count difference. (2) The log-scale interval is $-0.223\pm1.96(0.100)=[-0.419,-0.027]$; exponentiating gives approximately $[0.658,0.973]$. (3) $\exp(1.099+0.200)=\exp(1.299)\approx3.666$, or 3.67 events per clinic-period. (4) $\exp(1.099)$ is the mean at $b_j=0$; averaging over normally distributed random intercepts multiplies it by $\exp(\tau^2/2)$, giving $3.001\exp(0.080)\approx3.251$ events per clinic-period for the illustrative variance estimate.