Marginal Effects and Population-Averaged Estimates · Tutorial 403 of 698

The Integrated Mean In A GLMM

Express and evaluate a GLMM population-averaged response mean by integrating conditional means over the model’s random-effects distribution.

Advanced 12 min read

What You'll Learn

Define an integrated mean for a specified covariate pattern and random-effects distribution
Distinguish an integrated mean from a conditional mean at a particular random-effect value
Evaluate discrete random-effects integrals by weighting conditional means
Approximate continuous-distribution integrals using quadrature or simulation
Identify assumptions, numerical checks, and common interpretation errors

Why Integrate Over Random Effects?

A GLMM describes a response mean conditional on both observed predictors and random effects. That conditional mean answers a question about a specified random-effect value. A population-averaged mean asks a different question: what response mean does the model imply across the random-effects population specified by the model? To answer it, average the conditional mean over that distribution.

The integration is part of defining the estimand, not a change to the fitted conditional model. The covariate values and the distribution over which to average must be stated. Changing either can change the resulting mean. This tutorial focuses on the integral itself and its evaluation; the distinction between conditional and population-averaged estimands is the relevant earlier material.

Definition: Suppose the GLMM conditional mean at covariate values $x$ and random-effect value $b$ is $\mu(x,b)$, with $g\{\mu(x,b)\}=x^\mathsf{T}\beta+z^\mathsf{T}b$. For a specified random-effects distribution $F_b$, the integrated mean at those covariate values is $$ \bar{\mu}(x,z)=\int g^{-1}\!\left(x^\mathsf{T}\beta+z^\mathsf{T}b\right)\,dF_b(b). $$ It is the average of the conditional response mean over the specified distribution of $b$. If $F_b$ has density $f_b$, this is $\int g^{-1}(x^\mathsf{T}\beta+z^\mathsf{T}b)f_b(b)\,db$; for a discrete distribution, it is a weighted sum.

Model and Conditions for the Integral

For observation $j$ in cluster $i$, a common GLMM specification is $g\{E(Y_{ij}\mid b_i)\}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_i$. Here $x_{ij}$ and $z_{ij}$ are the fixed- and random-effect design vectors, $\beta$ is the fixed-effect coefficient vector, and $b_i$ is a random-effects vector with specified distribution $F_b(\,\cdot\,;\theta)$. The inverse link gives the conditional mean: $$ \mu_{ij}(b_i)=E(Y_{ij}\mid b_i) =g^{-1}\!\left(x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_i\right). $$

For an integrated mean, hold the covariate pattern fixed and average this conditional mean over $F_b$. The random-effects distribution must be a proper distribution, the inverse link must produce valid conditional means on its support, and the integral must exist. A sufficient integrability condition is $\int |\mu_{ij}(b)|\,dF_b(b)<\infty$. For bounded response means, such as binary-outcome probabilities, finiteness is automatic. For unbounded means it needs attention; the link and random-effects distribution can matter.

The target distribution is a substantive part of the definition. In a model with one random intercept, it may be the model’s distribution of cluster intercepts. With random slopes, it is the joint distribution of the intercept and slopes, including their dependence. The integral does not, by itself, specify an empirical weighting scheme over clusters or covariate patterns: those must be set separately if the target is an average over more than one covariate pattern.

Discrete Distributions: Weighted Conditional Means

When $b$ takes values $b_1,\ldots,b_K$ with probabilities $w_1,\ldots,w_K$, where $w_k\geq 0$ and $\sum_k w_k=1$, the integral is exactly $$ \bar{\mu}(x,z)=\sum_{k=1}^{K}w_k\,g^{-1}\!\left(x^\mathsf{T}\beta+z^\mathsf{T}b_k\right). $$ Each term is a conditional mean at one random-effect value; its weight is that value’s probability under the specified distribution. The result remains on the response-mean scale.

Worked Example: A Logistic Random Intercept With Three Support Points

Suppose a hypothetical binary-outcome model has logit link, fixed intercept $\beta_0=0.4$, and random intercept $b$ with probabilities $1/6$ at $-1$, $2/3$ at $0$, and $1/6$ at $1$. Find the integrated event probability for a covariate pattern with no additional fixed or random-slope contribution. The probabilities are nonnegative and sum to one, and the logistic inverse link returns values between zero and one, so the integral exists.

For the logit link, $g^{-1}(\eta)=1/(1+\exp(-\eta))$. The three linear predictors are $0.4-1=-0.6$, $0.4+0=0.4$, and $0.4+1=1.4$. Their conditional probabilities are approximately $0.3543$, $0.5987$, and $0.8022$, respectively. Therefore,

\begin{aligned} \bar{\mu} &=\frac{1}{6}\frac{1}{1+\exp(0.6)} +\frac{2}{3}\frac{1}{1+\exp(-0.4)} +\frac{1}{6}\frac{1}{1+\exp(-1.4)}\\ &\approx \frac{1}{6}(0.354344) +\frac{2}{3}(0.598688) +\frac{1}{6}(0.802184)\\ &\approx 0.059057+0.399125+0.133697\\ &\approx 0.5919. \end{aligned}

The model’s integrated event probability is about $0.5919$, or $59.19\%$, for this covariate pattern and this random-intercept distribution. The conditional probability at $b=0$ is about $0.5987$; it is not the integrated probability.

Continuous Distributions and Numerical Evaluation

When random effects have a continuous distribution, the same definition applies, but the integral may not have a closed-form solution. For a one-dimensional random intercept with density $f_b$, the target is $$ \bar{\mu}(x)=\int_{-\infty}^{\infty}g^{-1}(x^\mathsf{T}\beta+b)f_b(b)\,db. $$ A numerical method replaces the integral with a weighted sum or an average over draws. For example, quadrature uses selected evaluation points $b_k$ and weights $w_k$ chosen to approximate integration under the specified distribution: $$ \bar{\mu}(x)\approx\sum_{k=1}^{K}w_k\,g^{-1}(x^\mathsf{T}\beta+b_k). $$ For Monte Carlo integration, draw $b^{(1)},\ldots,b^{(M)}$ from the specified random-effects distribution and use $$ \bar{\mu}(x)\approx\frac{1}{M}\sum_{m=1}^{M}g^{-1}(x^\mathsf{T}\beta+b^{(m)}). $$

These approximations evaluate the model-based integral. Their random-effect values are drawn from the target distribution, rather than being restricted to estimated random effects for clusters already observed. With Monte Carlo integration, the simulation average also has numerical error; increasing $M$ and checking the stability of the estimate are practical safeguards.

Worked Example: A Two-Point Random Intercept Under a Complementary Log-Log Link

Consider a hypothetical model whose conditional event probability is $h(\eta)=1-\exp\{-\exp(\eta)\}$. At a fixed covariate pattern, let the linear predictor be $0.3+b$, where $b=0$ with probability $0.8$ and $b=-1$ with probability $0.2$. These probabilities define a proper distribution, and $h$ is between zero and one, so the integrated probability is well-defined.

At $b=0$, the conditional probability is $h(0.3)=1-\exp\{-\exp(0.3)\}\approx 0.7407$. At $b=-1$, it is $h(-0.7)=1-\exp\{-\exp(-0.7)\}\approx 0.3914$. Weight the conditional probabilities by their specified probabilities, using the unrounded values for the final calculation:

\begin{aligned} \bar{\mu} &=0.8\left[1-\exp\{-\exp(0.3)\}\right] +0.2\left[1-\exp\{-\exp(-0.7)\}\right]\\ &\approx 0.8(0.740723)+0.2(0.391395)\\ &\approx 0.592579+0.078279\\ &\approx 0.6709. \end{aligned}

The integrated event probability is approximately $0.6709$, or $67.09\%$. The arithmetic can be checked directly: $\exp(-0.496585)\approx0.6086$, so $h(-0.7)\approx1-0.6086=0.3914$. Keeping additional digits until the final rounding avoids a discrepancy caused by rounding intermediate terms.

Joint Random Effects and a Covariate-Specific Mean

When $b$ is a vector, the integral is over its joint distribution. The covariate pattern determines how each component enters the linear predictor through $z$. It is generally not enough to integrate each random-effect component separately while ignoring dependence: the required distribution is the joint $F_b$ specified by the model.

Worked Example: Integrating a Random Intercept and Slope

For a hypothetical logistic model at $x=1$, let the linear predictor be $-0.4+0.7x+b_0+b_1x$. Suppose the joint random-effects distribution has three support points: $(b_0,b_1)=(0,0)$ with probability $0.5$, $(-0.2,-0.3)$ with probability $0.25$, and $(0.2,0.3)$ with probability $0.25$. The weights sum to one, and the logistic inverse link gives finite probabilities at all three points.

At $x=1$, the fixed part is $-0.4+0.7=0.3$. The random-effect contributions $b_0+b_1x$ are $0$, $-0.5$, and $0.5$, giving linear predictors $0.3$, $-0.2$, and $0.8$. The corresponding logistic means are approximately $0.574443$, $0.450166$, and $0.689974$. Thus,

\begin{aligned} \bar{\mu}(1) &=0.5\frac{1}{1+\exp(-0.3)} +0.25\frac{1}{1+\exp(0.2)} +0.25\frac{1}{1+\exp(-0.8)}\\ &\approx 0.5(0.574443)+0.25(0.450166)+0.25(0.689974)\\ &\approx 0.287221+0.112542+0.172494\\ &\approx 0.5723. \end{aligned}

For this specified joint distribution and $x=1$, the integrated event probability is about $0.5723$, or $57.23\%$. The same calculation at a different $x$ would use different random-effect contributions and can give a different integrated mean.

Evaluation, Diagnostics, and Alternatives

A useful evaluation sequence is to state the target, identify the distribution being integrated over, compute the conditional mean at each required random-effect value, and then take the weighted integral. For a continuous distribution, replace the exact integral by a suitable numerical approximation and check that it is accurate enough for the intended use.

Specify the target covariates.
Record the fixed-effect pattern $x$ and the random-effect design $z$. State any covariate values held fixed.

Specify the random-effects distribution.
Use the intended $F_b$, including its joint structure when there is more than one random effect.

Evaluate and average conditional means.
Apply $g^{-1}$ to the linear predictor at each integration point, then integrate with the appropriate probabilities, density, or numerical weights.

Check the result.
Confirm weights and probabilities are valid, the mean is in the response family’s allowable range, and numerical approximations are stable under increased integration accuracy.

For quadrature, a practical diagnostic is to recompute the mean with more integration points and compare the results at the precision needed for reporting. For Monte Carlo integration, repeat the calculation with a larger number of draws or independent batches; substantial variation indicates that more draws may be needed. Also check that the implemented random-effects distribution matches the model specification, particularly its scale and covariance.

The main alternatives answer different questions. A conditional prediction at $b=0$ evaluates the mean for a zero random effect; a cluster-specific prediction using an estimated random effect targets that cluster’s fitted conditional mean. Neither is the same operation as integrating over the population distribution. A marginal model fitted directly to population-averaged responses is another approach, with its own modeling and inferential framework. The GLMM integral instead derives the population-averaged mean from the specified conditional model and random-effects distribution.

Common Mistakes and Limitations

Substituting a typical random effect for integration. Setting $b=0$, or using a single estimated cluster effect, gives one conditional mean, not an average over $F_b$.
Averaging linear predictors instead of response means. The defined integral averages $g^{-1}(x^\mathsf{T}\beta+z^\mathsf{T}b)$. Averaging the linear predictor and applying the inverse link is a different calculation.
Using the wrong random-effects distribution. The result depends on the specified distribution and, with multiple effects, their joint dependence. A convenient substitute distribution changes the target calculation.
Ignoring numerical error. A plausible-looking numerical value is not proof of accurate integration. Check convergence or simulation variability, and avoid rounding intermediate conditional means too early.

The integrated mean inherits the limitations of the fitted conditional model and the assumed random-effects distribution. It is not automatically a universally representative population quantity: its interpretation is tied to the covariate pattern and the distribution used for averaging. Estimating $\beta$ and the distribution parameters also introduces statistical uncertainty beyond numerical integration error; an integration routine’s accuracy does not quantify that inferential uncertainty.

Key takeaway: At fixed covariates, a GLMM integrated mean is the random-effects-distribution average of the inverse-linked conditional linear predictor. State the target distribution explicitly, integrate response means rather than linear predictors, and separate numerical accuracy from uncertainty in the fitted model.

Check Your Understanding

Answer each question, then compare with the solutions below.

A random intercept takes values $-1$ and $1$, each with probability $0.5$. With a logit-link intercept of $0$, write the integrated mean as a weighted sum.
In a random-intercept model, why is the conditional mean at $b=0$ not generally the integrated mean?
For a random intercept and slope, which distribution is required in the integral?
A Monte Carlo estimate changes noticeably when the number of draws is doubled. What should be checked or changed?

Solutions. (1) The linear predictors are $-1$ and $1$, so $\bar{\mu}=0.5\{1/(1+\exp(1))\}+0.5\{1/(1+\exp(-1))\}=0.5(0.2689)+0.5(0.7311)=0.5000$. (2) It evaluates one random-effect value, whereas the integrated mean averages conditional means over the full specified distribution. (3) The joint distribution of the intercept and slope, including their dependence. (4) Increase the number of draws, assess variability across independent batches, and verify that draws come from the intended distribution; the change signals that simulation error may still matter.

Back to Tutorials