Tutorials › Generalized Linear Mixed Models (GLMM) › Predicting Group Effects in a Poisson GLMM

Prediction of Random Effects · Tutorial 457 of 698

Predicting Group Effects in a Poisson GLMM

Learn to obtain and interpret a group’s conditional-mode random intercept in a Poisson GLMM, accounting for both its observed count and its exposure.

Advanced 10 min read

What You'll Learn

  • Write the Poisson random-intercept model with an exposure offset.
  • Derive the score equation whose solution is the group conditional mode.
  • Check why the score has a unique root and how to compute it numerically.
  • Interpret predictions for groups with high counts and with zero counts.
  • Explain how exposure and the random-effect variance affect shrinkage.
  • Diagnose calculation errors and distinguish group-specific from population-averaged predictions.

Why a Group Count Does Not Directly Determine Its Effect

In a Poisson generalized linear mixed model, a group’s observed count reflects both its underlying group effect and the amount of exposure it had. A group with twice as many events might simply have twice the follow-up time. Even after accounting for exposure, a raw group-specific rate can be unstable when the count is small. Conditional-mode prediction combines the group’s data with the model’s random-effects distribution to produce a group-specific estimate.

The target here is the conditional mode of one group’s random intercept, given that group’s count data and fixed model parameters. It is a prediction on the log-rate scale, not a population-average effect and not a direct estimate of a probability.

Definition: For groups \(j=1,\ldots,J\), let \(Y_{ij}\) be the count for observation \(i\) in group \(j\), and let \(e_{ij}>0\) be its exposure. A Poisson log-link random-intercept model is $$ Y_{ij}\mid b_j \sim \operatorname{Poisson}(\mu_{ij}), \qquad \log(\mu_{ij})=\log(e_{ij})+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_j, \qquad b_j\sim N(0,\tau^2). $$ Conditional on the random effects, the counts are independent. The group effect \(b_j\) is independent of the covariates, and \(\tau^2>0\). The conditional-mode prediction for group \(j\) is the value of \(b_j\) maximizing its conditional density given that group’s observed counts, with \(\boldsymbol{\beta}\) and \(\tau^2\) held fixed.

The exposure term is an offset: increasing exposure multiplies the expected count while leaving the underlying rate structure unchanged. For a given group, define the total observed count and its baseline expected count, before applying the group effect, by \(y_{+j}=\sum_i y_{ij}\) and \(A_j=\sum_i e_{ij}\exp(\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta})\). The value \(A_j\) accounts for both exposures and any differences in fixed covariates. It is not generally just the number of observations.

Deriving the Conditional Mode

Conditional on a proposed value \(b_j\), the contribution of group \(j\)’s counts to the log likelihood, apart from terms not involving \(b_j\), is \(y_{+j}b_j-A_j\exp(b_j)\). The normal random-effects density contributes the penalty \(-b_j^2/(2\tau^2)\). Thus the log of the conditional density, up to an additive constant, is

$$ \ell_j(b_j)=y_{+j}b_j-A_j\exp(b_j)-\frac{b_j^2}{2\tau^2}. $$

Differentiating gives the score equation. Its root, rather than the value that would exactly match the observed rate, is the conditional mode.

$$ U_j(b_j)=\frac{d\ell_j}{db_j} =y_{+j}-A_j\exp(b_j)-\frac{b_j}{\tau^2}=0. $$

There is exactly one solution. The second derivative is \(\ell_j''(b_j)=-A_j\exp(b_j)-1/\tau^2<0\), so the log conditional density is strictly concave. Also, \(U_j(b_j)\) tends to positive infinity as \(b_j\) tends to negative infinity and to negative infinity as \(b_j\) tends to positive infinity. Since \(U_j\) is continuous and strictly decreasing, it crosses zero exactly once. That root is the unique maximizer.

In general the root is found numerically. Newton’s method, starting at a value \(b_j\), uses

$$ b_j^{\mathrm{new}} =b_j+\frac{y_{+j}-A_j\exp(b_j)-b_j/\tau^2} {A_j\exp(b_j)+1/\tau^2}. $$

The denominator is positive, so the update is well-defined; the updated value need not be positive. Convergence should be checked by evaluating the score at the final value, not merely by noting that an algorithm stopped.

Worked Example: A Group With a High Count

Suppose a hypothetical group has total count \(y_{+j}=28\), baseline expected count \(A_j=20\), and random-intercept variance \(\tau^2=2\). Find its conditional mode. Treat these model quantities as fixed for the calculation.

The score equation is \(28-20\exp(b_j)-b_j/2=0\). Solving numerically gives \(b_j\approx0.33055\). Check the root by substitution: \(20\exp(0.33055)\approx27.83467\) and \(0.33055/2=0.165275\). Therefore \(28-27.83467-0.165275\approx0.00006\).

The distinction between \(\tau^2=2\) and the penalty coefficient matters. For the commonly used setting \(1/\tau^2=2\), equivalently \(\tau^2=0.5\), the score equation is \(28-20\exp(b_j)-2b_j=0\). Its root is about \(0.31380\): \(20\exp(0.31380)\approx27.3723\), \(2(0.31380)=0.6276\), and the score is \(28-27.3723-0.6276\approx0.0001\). For this setting, \(\exp(b_j)\approx1.3686\), so the predicted group rate is about 1.369 times its baseline rate.

This calculation also illustrates why stating the variance precisely is essential. The conditional mode depends on \(\tau^2\), not just on the observed count and baseline expectation.

For clarity in the next examples, all calculations state the variance and use the corresponding penalty coefficient \(1/\tau^2\) explicitly. When the baseline covariates and exposures are equal within a group, \(A_j\) is its exposure multiplied by its baseline rate. The predicted conditional mean count is then \(A_j\exp(\hat b_j)\).

Worked Example: A Group With No Events

A hypothetical group has \(y_{+j}=0\), \(A_j=5\), and \(\tau^2=0.8\). Its score equation is \(0-5\exp(b_j)-b_j/0.8=0\), or equivalently \(5\exp(b_j)+1.25b_j=0\).

Numerical solution gives \(\hat b_j\approx-1.2022\). At this value, \(5\exp(-1.2022)\approx1.5027\) and \(1.25(-1.2022)=-1.5028\), so the score is approximately zero to the displayed precision. The conditional predicted mean count is \(5\exp(-1.2022)\approx1.503\), not zero.

The negative effect reflects the absence of events relative to the baseline expectation of 5. It does not drive the predicted mean all the way to zero because the normal random-effects density penalizes large departures from zero. For a group with zero events, the mode is negative, but it remains finite.

Worked Example: More Exposure Changes the Comparison

Suppose a hypothetical group has baseline rate \(0.4\) events per exposure unit, exposure \(20\) units, observed count \(12\), and \(\tau^2=0.5\). The baseline expected count is \(A_j=20(0.4)=8\). The score equation is \(12-8\exp(b_j)-2b_j=0\).

Numerical solution gives \(\hat b_j\approx0.34606\). Substitution gives \(8\exp(0.34606)\approx11.3079\) and \(2(0.34606)=0.69212\), so the score is approximately \(12-11.3079-0.69212\approx0.0000\). The conditional predicted mean count is about 11.308.

The observed rate is \(12/20=0.600\) events per exposure unit, 1.5 times the baseline rate. Without the random-effects penalty, matching that rate would give \(\log(0.600/0.400)=\log(1.5)\approx0.4055\). The conditional mode is smaller, about \(0.3461\), and the predicted rate is \(0.4\exp(0.34606)\approx0.5654\) events per exposure unit. Thus the estimate accounts for the exposure and is pulled toward the population’s zero-centered random-intercept distribution.

How Counts, Exposure, and Variance Shape the Prediction

The score equation makes the roles of the data and model visible. A larger count tends to move the root upward; a larger baseline expected count \(A_j\), with the observed count held fixed, tends to move it downward. Exposure matters through \(A_j\): a count of 12 is more evidence of a high rate over 20 exposure units than over 200 units, all else equal. The variance controls the strength of the normal penalty. A smaller \(\tau^2\) makes \(1/\tau^2\) larger and pulls the mode more strongly toward zero; a larger \(\tau^2\) weakens that penalty.

When \(y_{+j}=A_j\), \(b_j=0\) solves the score equation. When \(y_{+j}>A_j\), the score at zero is positive, so the unique root is positive; when \(y_{+j}<A_j\), it is negative. These sign statements follow from \(U_j(0)=y_{+j}-A_j\) and the fact that \(U_j\) decreases strictly.

For a practical computation, first calculate \(y_{+j}\), \(A_j\), and \(1/\tau^2\) using the model’s exposure and fixed predictors. Then solve the scalar score equation and verify its residual. Finally, obtain fitted conditional means using \(\hat\mu_{ij}=e_{ij}\exp(\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+\hat b_j)\). Each mean is a count for that observation’s exposure; the rate per exposure unit is \(\exp(\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+\hat b_j)\).

Diagnostics, Alternatives, and Limitations

A useful numerical diagnostic is to evaluate \(U_j(\hat b_j)\). A substantial nonzero residual means the reported value does not solve the mode equation to the stated precision. Also check that the calculated \(A_j\) includes every exposure and fixed covariate contribution. Using the number of records instead of exposure, or omitting a log-exposure offset, changes the model being predicted. If a calculation produces a positive mode when \(y_{+j}<A_j\), or a negative mode when \(y_{+j}>A_j\), check the score, the variance, and the exposure construction.

A raw rate comparison would set a group effect to \(\log(y_{+j}/A_j)\) when \(y_{+j}>0\). That value matches the group’s observed rate exactly but ignores the random-effects distribution; it is undefined when the count is zero. The conditional mode instead balances fit to the group data against the normal penalty. A population-averaged prediction answers a different question: it averages over the distribution of random effects rather than conditioning on a particular group’s data. With a log link and normal random intercept, integrating over \(b_j\) changes the mean rate by the factor \(\exp(\tau^2/2)\), so the population average is not simply the prediction obtained by setting \(b_j=0\).

The conditional mode is a point prediction, not an interval and not a statement that the group’s true effect equals the predicted value. Its calculation conditions on the specified model parameters. It also relies on the Poisson mean-variance relationship, the log link, the exposure specification, and the normal random-intercept distribution. If the data show substantial overdispersion beyond the model, the exposure is recorded incorrectly, or the random-effects structure is unsuitable, the numerical root can still be calculated but may not represent a useful group prediction.

Common pitfall: Treating the observed group rate as the predicted group effect. The raw log rate ratio \(\log(y_{+j}/A_j)\) ignores the random-effects distribution and fails at zero counts. The conditional mode instead solves \(y_{+j}-A_j\exp(b_j)-b_j/\tau^2=0\), so it uses both the group data and the model’s penalty.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. If \(y_{+j}=A_j\), what is the conditional mode?
  2. For \(y_{+j}=0\), \(A_j=4\), and \(\tau^2=1\), state the sign of the mode and explain why it is not negative infinity.
  3. A group has count 10 and baseline expected count 8. Is its conditional mode positive, zero, or negative?
  4. Why should two groups with the same count but different exposures not automatically receive the same prediction?

Solutions. (1) The mode is zero, since substituting \(b_j=0\) gives score \(y_{+j}-A_j=0\), and the root is unique. (2) The mode is negative because \(U_j(0)=0-4<0\) and the score decreases; it is finite because the score tends to positive infinity as \(b_j\) tends to negative infinity and crosses zero once. (3) It is positive because \(U_j(0)=10-8=2>0\), and the decreasing score reaches zero at a positive value. (4) Exposure contributes to \(A_j\), the baseline expected count. The same count relative to a larger exposure is weaker evidence for a high group rate, so the score equation—and generally its root—differs.

Clinical Biostats