Pitfall: Averaging Linear Predictors Instead Of Responses

Why the Scale of Averaging Matters

In a generalized linear mixed model, the linear predictor is on a scale transformed by the link function. For a binary outcome with a logit link, a linear predictor of $0.8$ corresponds to a probability of about $0.690$, not to a probability of $0.8$. When random effects vary across clusters, the probability also varies across clusters. A population-averaged probability must therefore average the probabilities represented by those different linear predictors.

A tempting shortcut is to average the linear predictors first and then apply the inverse link. This can produce a valid probability, but it generally describes a different quantity: the inverse-link transformation of the average predictor. Nonlinearity means that this is usually not the average response probability. The distinction matters whenever results are described as population-averaged predictions or contrasts.

Model and Target Quantity

Consider a binary outcome $Y_{ij}$ for individual $i$ in cluster $j$. For a random-intercept logistic GLMM, suppose

Y_{ij}\mid b_j \sim \operatorname{Bernoulli}(p_{ij}(b_j)), \qquad \operatorname{logit}\{p_{ij}(b_j)\}=\eta_{ij}+b_j, \qquad b_j\sim N(0,\sigma_b^2).

Here $\eta_{ij}=x_{ij}^{\mathsf T}\beta$ is the fixed-effects linear predictor at the specified covariates, $b_j$ is the cluster-specific random intercept, and $\operatorname{logit}(p)=\log\{p/(1-p)\}$. The inverse-logit function is $\operatorname{expit}(u)=1/(1+\exp(-u))$, so the conditional probability is $\operatorname{expit}(\eta_{ij}+b_j)$. The model assumes that, conditional on the random effect and covariates, the outcome is Bernoulli with this probability, and that the random intercepts follow the stated normal distribution. When using the usual random-intercept specification, the random effects are independent across clusters and independent of the included covariates.

Definition: At fixed covariates $x$, the random-effects population-averaged probability is $$ \bar p(x)=E_{b}\left[\operatorname{expit}\{x^{\mathsf T}\beta+b\}\right] =\int \operatorname{expit}(x^{\mathsf T}\beta+b)f_b(b)\,db. $$ It averages the response probability over the specified random-effects distribution. In general, it is not $\operatorname{expit}\{x^{\mathsf T}\beta+E(b)\}$.

Because $E(b)=0$ in this model, $\operatorname{expit}(x^{\mathsf T}\beta)$ is the probability for a cluster with random intercept $b=0$. It is not, in general, the probability averaged across the random-effects population. With a random slope, the same principle applies: first form the conditional probability for each random-effects value, then average those probabilities under the relevant joint random-effects distribution.

Why the Shortcut Fails

Let $g(u)=\operatorname{expit}(u)$, and write $\eta=x^{\mathsf T}\beta$. Averaging the predictors and applying the inverse link gives $g\{E(\eta+b)\}=g(\eta)$, since $E(b)=0$. The desired population average is $E\{g(\eta+b)\}$. These quantities would be equal for every distribution of $b$ only if the inverse link were linear over the relevant range. The logistic function is nonlinear, so equality is not generally expected.

A second-order Taylor expansion around $\eta$ makes the discrepancy visible. For a mean-zero random effect with variance $\sigma_b^2$, the first-order term averages to zero, giving the approximation

E\{g(\eta+b)\} \approx g(\eta)+\frac{1}{2}\sigma_b^2g''(\eta), \qquad g''(\eta)=g(\eta)\{1-g(\eta)\}\{1-2g(\eta)\}.

Thus the approximation difference depends on both the random-effects variance and the location of $\eta$. Below probability $0.5$, the logistic curve is locally convex; above $0.5$, it is locally concave. There is no universal rule that the average probability must be larger or smaller than the inverse link of the average predictor. The Taylor expression is a justification for the discrepancy, not a substitute for integration when an accurate marginal probability is needed.

Calculating the Response-Scale Average

For a normal random intercept, the required integral is a logistic-normal integral. It usually has no elementary closed form, so it can be evaluated using numerical quadrature or simulation. A transparent approximation uses three-point Gaussian quadrature for a standard normal variable: use nodes $-\sqrt{3},0,\sqrt{3}$ with weights $1/6,2/3,1/6$. For $b\sim N(0,\sigma_b^2)$, multiply the nodes by $\sigma_b$. This gives

\bar p(x)\approx \frac{1}{6}g(\eta-\sqrt{3}\sigma_b) +\frac{2}{3}g(\eta) +\frac{1}{6}g(\eta+\sqrt{3}\sigma_b).

This is an approximation to the normal integral, not a claim that a normal random effect takes only three possible values. More quadrature points or simulation can be used when greater numerical accuracy is needed.

Worked Example: One Covariate Pattern

Suppose a hypothetical binary random-intercept GLMM has fixed-effects predictor $\eta=0.8$ at a specified covariate pattern and random-intercept standard deviation $\sigma_b=1.0$. Compare the zero-random-effect probability with the approximate population-averaged probability.

The zero-random-effect prediction is $g(0.8)=1/\{1+\exp(-0.8)\}=0.6900$. For quadrature, the nonzero nodes are $\pm\sqrt{3}(1.0)=\pm1.732$. The corresponding probabilities are $g(0.8-1.732)=g(-0.932)=0.2825$ and $g(0.8+1.732)=g(2.532)=0.9264$. The center-node probability remains $0.6900$.

Substituting into the three-point formula gives

\bar p(x)\approx \frac{1}{6}(0.2825)+\frac{2}{3}(0.6900)+\frac{1}{6}(0.9264) =0.04708+0.4600+0.1544 =0.6615.

The approximate population-averaged probability is $0.6615$, or about 66.2%, whereas $0.6900$ is the probability for a cluster with $b=0$. Applying the inverse link to the mean predictor returns the latter, not the former. Both are probabilities, but they answer different questions.

Worked Example: A Binary Exposure Contrast

Suppose a hypothetical model has control predictor $\eta_0=-0.8$, exposure predictor $\eta_1=0.2$, and shared random-intercept standard deviation $\sigma_b=0.8$. Compute the conditional and population-averaged probability contrasts at these covariates. Assume the model and normal random-intercept assumptions stated above.

At $b=0$, the control probability is $g(-0.8)=0.3100$, and the exposure probability is $g(0.2)=0.5498$. The conditional probability difference is therefore $0.5498-0.3100=0.2398$, or 24.0 percentage points.

For the marginal probabilities, the quadrature nodes are $0$ and $\pm\sqrt{3}(0.8)=\pm1.3856$. For control, the probabilities at the negative, center, and positive nodes are $g(-2.1856)=0.1011$, $g(-0.8)=0.3100$, and $g(0.5856)=0.6424$. Hence

\bar p_0\approx \frac{1}{6}(0.1011)+\frac{2}{3}(0.3100)+\frac{1}{6}(0.6424) =0.3306.

For exposure, the corresponding linear predictors are $-1.1856$, $0.2$, and $1.5856$, giving probabilities $0.2340$, $0.5498$, and $0.8300$. Thus

\bar p_1\approx \frac{1}{6}(0.2340)+\frac{2}{3}(0.5498)+\frac{1}{6}(0.8300) =0.5439.

The population-averaged probability difference is approximately $0.5439-0.3306=0.2133$, or 21.3 percentage points. It is not the conditional difference of 24.0 percentage points. The calculation averages each condition's response probabilities over the same random-intercept distribution before taking their difference.

Worked Example: Averaging Across a Finite Set of Clusters

Suppose a hypothetical prediction exercise gives four equally weighted clusters random-intercept values $-1.2,-0.4,0.3,1.3$, whose mean is zero. At a fixed covariate pattern, let the fixed predictor be $\eta=0.5$. Find the average predicted probability across these four clusters and compare it with applying the inverse link to the average predictor.

The four cluster-specific predictors are $-0.7,0.1,0.8,1.8$, with inverse-logit probabilities $0.3318,0.5250,0.6900,0.8581$. Their equally weighted response average is

\frac{0.3318+0.5250+0.6900+0.8581}{4} =\frac{2.4049}{4}=0.6012.

The average predictor is $\{(-0.7)+0.1+0.8+1.8\}/4=0.5$, so its inverse link is $g(0.5)=0.6225$. It differs from the finite-cluster response average, $0.6012$. This finite average describes these four equally weighted cluster predictions; it need not equal the normal random-effects population integral. The example isolates the scale error without treating a finite set of cluster effects as if it were the assumed normal distribution.

Choosing the Correct Average and Checking the Calculation

Before calculating, specify what population the average represents. A normal random-effects integral averages over the model's random-effects distribution at fixed covariates. An equally weighted average over observed clusters gives each listed cluster the same weight. An individual-weighted average gives clusters with more individuals greater influence. These targets need not agree. The response-scale correction addresses the order of averaging and transformation; it does not by itself decide which covariate distribution or cluster weights are appropriate.

A practical check is to write down the intended estimand in words and as an equation before producing a number. For a population-averaged probability, the inverse link belongs inside the expectation:

E_b\left[g\{\eta(x)+b\}\right], \qquad\text{not}\qquad g\left(E_b[\eta(x)+b]\right).

For a contrast between two covariate conditions, calculate each average probability using the same target random-effects distribution, then subtract them. If covariates are also being averaged over a target population, that averaging must also be defined and applied to response probabilities. A clear label such as “probability at $b=0$” prevents a conditional prediction from being reported as though it were population-averaged.

Numerical diagnostics include checking that the quadrature weights sum to one, that each calculated probability lies between zero and one, and that increasing the number of quadrature points leaves the answer essentially unchanged. A large change as quadrature is refined indicates that the initial approximation may be inadequate. Simulation offers another check: draw many random effects from the specified distribution, calculate the inverse-link probability for each draw, and average those probabilities. The simulation average should approach the integral as the number of draws increases.

Common Mistake, Alternatives, and Limits

Wrong reasoning: “The random intercept has mean zero, so averaging it out leaves the fixed predictor; therefore the inverse link of the fixed predictor is the population-averaged probability.” The first part is correct for the mean predictor, but the conclusion moves the inverse link outside a nonlinear expectation.
Correction: For a population average on the response scale, calculate $g(\eta+b)$ for random-effects values and average those probabilities under the target distribution.
Alternative target: If the question concerns a cluster with $b=0$, then $g(\eta)$ is appropriate. It should be described as a zero-random-effect or reference-cluster prediction, not as an average over clusters.
Alternative calculation: Numerical quadrature and simulation both approximate the response-scale integral. They are computational alternatives for the same target when they use the same random-effects distribution and covariates.

The response-scale integral depends on the fitted fixed effects, the random-effects distribution, and its variance. If that distribution is misspecified, the resulting marginal probability may not represent the intended population. Quadrature accuracy does not correct model misspecification, and simulation error does not disappear unless the simulation is sufficiently large. Finally, a contrast between two marginal probabilities is not automatically a causal effect: its interpretation also depends on the study design and assumptions supporting the comparison.

Key takeaway: Transforming an average predictor gives a prediction at the average predictor. A population-averaged GLMM probability instead averages inverse-link predictions over the random-effects distribution. State which target is intended before choosing the calculation.

Check Your Understanding

Answer each question, then compare with the brief solutions.

In a logistic random-intercept GLMM with $E(b)=0$, what does $g(\eta)$ represent, and why is it not generally the population-averaged probability?
For $\eta=0$ and $\sigma_b=0$, what are the zero-random-effect and population-averaged probabilities?
Why should probabilities be averaged before calculating a difference between two population-averaged conditions?
A researcher averages cluster-specific predictors and then applies the inverse logit. Name one additional check needed even after correcting the order of operations.

Solutions. (1) $g(\eta)$ is the probability at $b=0$; the population average is $E_b[g(\eta+b)]$, which generally differs because $g$ is nonlinear. (2) Both equal $g(0)=0.5000$, since $b=0$ with probability one when its variance is zero. (3) Each condition has its own response probability at each random-effects value, so average each condition's probabilities over the same target distribution before subtracting. (4) Specify the target population and weighting, verify the random-effects distribution and covariates match that target, or check numerical integration accuracy by refining quadrature or comparing with simulation.