From Cluster-Specific Lines to a Population Curve
When a covariate changes, a GLMM with random slopes allows clusters to have different response patterns: their curves can start at different levels and change at different rates. A single prediction for a cluster with random effects set to zero describes a reference cluster, not the average response across a population of clusters. A population-averaged curve instead averages the conditional response over the random-effects distribution at each covariate value.
The curve is useful when the target is the expected response for a randomly selected cluster from the model’s target population, at specified covariate values. Its shape depends on the inverse link and on how random effects enter the linear predictor. In general, it is not obtained by averaging the random effects first and then applying the inverse link.
For a random intercept and random slope, write the conditional model as $$ g\{\mu(x,b)\}=\beta_0+\beta_1x+b_0+b_1x, $$ where \(g\) is the link, \(\beta_0\) and \(\beta_1\) are fixed effects, and \(b=(b_0,b_1)^\mathsf{T}\) represents cluster-specific deviations. The standard model assumes \(b\) follows a specified distribution, commonly a bivariate normal distribution with mean zero and covariance matrix \(G\). Conditional on \(b\) and the covariates, responses follow the chosen GLM response distribution and are independent under the model. The marginal curve is therefore $$ m(x)=\int g^{-1}\{\beta_0+\beta_1x+b_0+b_1x\}\,dF_b(b). $$
This definition treats the random effects as varying across clusters in the target population. It does not average over a distribution of covariates: \(x\) is fixed separately at each point on the curve. If other covariates are present, they must also be specified or averaged over using a separately defined target distribution.
Constructing the Curve
A direct construction evaluates the conditional mean at a chosen value of \(x\) for many random-effects values, then averages those means using the model’s random-effects distribution. Repeating this at a grid of covariate values yields the curve. Numerical integration, quadrature, or Monte Carlo draws can approximate the integral. The averaging must be performed after applying the inverse link to each cluster-specific linear predictor.
Specify the values of \(x\) at which the population response is to be reported, and fix any other predictors as required.
For each random-effects value \(b\), calculate \(g^{-1}(\beta_0+\beta_1x+b_0+b_1x)\) at each chosen \(x\).
Integrate or take a suitably weighted numerical average of the conditional means over \(F_b\), separately for each \(x\).
Give response units, covariate units, the target population, and whether the reported values are estimates or uncertainty intervals.
A short justification explains why the order matters. Let \(h=g^{-1}\), and define the cluster-specific linear predictor \(\eta(x,b)=\beta_0+\beta_1x+b_0+b_1x\). The law of total expectation gives $$ E(Y\mid x)=E_b\{E(Y\mid x,b)\}=E_b\{h(\eta(x,b))\}. $$ In general, this is not \(h(E_b\{\eta(x,b)\})\), since a nonlinear function does not generally commute with averaging. Since the random effects have mean zero, \(E_b\{\eta(x,b)\}=\beta_0+\beta_1x\), but applying the inverse link to that fixed-effects predictor gives only the conditional mean at \(b=0\), not the population average.
Worked Example: A Logistic Curve From Numerical Averaging
Suppose a hypothetical binary-outcome GLMM has conditional probability $$ P(Y=1\mid x,b)=\operatorname{expit}(0.4-0.6x+b_0+b_1x), \qquad \operatorname{expit}(u)=\frac{1}{1+e^{-u}}. $$ To make the calculation transparent, use four equally weighted random-effects values as a numerical approximation: \((-0.5,-0.4)\), \((-0.2,-0.1)\), \((0.2,0.1)\), and \((0.5,0.4)\). These are illustrative integration points, not observed clusters or a claim that the random-effects distribution has only four possible values. The same calculation applies to quadrature points or random draws with their appropriate weights.
At \(x=1\), the random-effect contributions \(b_0+b_1\) are \(-0.9,-0.3,0.3,0.9\). The fixed part is \(0.4-0.6=-0.2\), so the four linear predictors are \(-1.1,-0.5,0.1,0.7\). Their response probabilities are approximately \(0.2497,0.3775,0.5250,0.6682\). Thus the averaged probability is $$ \widehat m(1)=\frac{0.2497+0.3775+0.5250+0.6682}{4}=0.4551. $$
At \(x=0\), the random contributions are just the intercept deviations. The predictors are \(-0.1,0.2,0.6,0.9\); their probabilities are \(0.4750,0.5498,0.6457,0.7109\). Hence \(\widehat m(0)=2.3814/4=0.5954\). At \(x=2\), the random contributions \(b_0+2b_1\) are \(-1.3,-0.4,0.4,1.3\). The fixed part is \(-0.8\), giving predictors \(-2.1,-1.2,-0.4,0.5\), and probabilities \(0.1091,0.2315,0.4013,0.6225\). Therefore \(\widehat m(2)=1.3644/4=0.3411\).
The estimated curve decreases across these three covariate values: the average response probabilities are about \(0.5954\), \(0.4551\), and \(0.3411\) at \(x=0,1,2\), respectively. At \(x=1\), applying the inverse link to the average predictor instead would give \(\operatorname{expit}(-0.2)=0.4502\), not \(0.4551\). The difference illustrates why averaging conditional probabilities, rather than predictors, defines the population-averaged response.
When Integration Has a Closed Form
Some links and random-effects distributions yield an analytic curve. For a probit model, suppose \(Y\) is binary and $$ P(Y=1\mid x,b)=\Phi(\beta_0+\beta_1x+b_0+b_1x), $$ where \(\Phi\) is the standard normal cumulative distribution function and \(b\) is bivariate normal with mean zero and covariance matrix \(G\). Let \(v(x)\) be the variance of \(b_0+b_1x\). Since a standard normal variable plus an independent normal random effect is normal, averaging the conditional probit probability gives $$ m(x)=\Phi\left(\frac{\beta_0+\beta_1x}{\sqrt{1+v(x)}}\right). $$ This identity follows by introducing an independent standard normal variable \(Z\): the conditional probability is \(P\{Z\leq \beta_0+\beta_1x+b_0+b_1x\}\). Averaging over \(b\) gives the probability that \(Z-(b_0+b_1x)\leq\beta_0+\beta_1x\). The left side is normal with mean zero and variance \(1+v(x)\), which produces the formula.
Worked Example: A Probit Population-Averaged Curve
Consider a hypothetical probit model with \(\beta_0=-0.4\), \(\beta_1=0.7\), and independent normal random effects with variances \(0.36\) for the intercept and \(0.16\) for the slope. Thus \(v(x)=0.36+0.16x^2\). The marginal probability is $$ m(x)=\Phi\left(\frac{-0.4+0.7x}{\sqrt{1.36+0.16x^2}}\right). $$
At \(x=-1\), the numerator is \(-1.1\), and the denominator is \(\sqrt{1.52}=1.233\). Therefore \(m(-1)=\Phi(-0.8922)\approx0.1862\). At \(x=0\), the ratio is \(-0.4/\sqrt{1.36}=-0.3430\), giving \(m(0)\approx0.3658\). At \(x=1\), the ratio is \(0.3/\sqrt{1.52}=0.2433\), giving \(m(1)\approx0.5961\). These are population-averaged probabilities at the stated covariate values; they are not probabilities for a cluster whose random effects have been observed.
A second useful closed-form case is a log-link model with normal random effects. If the conditional mean is \(\exp(\beta_0+\beta_1x+b_0+b_1x)\), then the random contribution is normal with mean zero and variance \(v(x)\). The moment-generating function of a normal variable gives $$ m(x)=\exp\{\beta_0+\beta_1x+\tfrac{1}{2}v(x)\}. $$ This is a mean on the response scale, such as an expected count or rate, and requires the stated normal random-effects assumption. It is not obtained by merely exponentiating the fixed-effects predictor.
Worked Example: A Log-Link Mean Curve
Suppose a hypothetical Poisson log-link model has unit exposure, fixed effects \(\beta_0=0.2\) and \(\beta_1=0.3\), and independent normal random intercept and slope with variances \(0.25\) and \(0.04\). Then \(v(x)=0.25+0.04x^2\), so $$ m(x)=\exp\{0.2+0.3x+\tfrac12(0.25+0.04x^2)\}. $$
At \(x=0\), the exponent is \(0.2+0.125=0.325\), giving \(m(0)=e^{0.325}\approx1.384\) expected events per unit exposure. At \(x=1\), the exponent is \(0.2+0.3+0.145=0.645\), giving \(m(1)\approx1.906\). At \(x=2\), it is \(0.2+0.6+0.205=1.005\), giving \(m(2)\approx2.732\). The marginal curve is not simply the exponential of the fixed-effects line because the random-effect variance contributes to the integrated mean.
Uncertainty, Diagnostics, and Interpretation
In practice, the fixed effects and random-effects distribution are estimated. A curve calculated by inserting those estimates into the integral is a plug-in estimate of \(m(x)\). Uncertainty in the curve reflects uncertainty in both the fixed effects and the random-effects distribution; numerical integration error is a separate computational issue. Pointwise confidence intervals describe uncertainty at each specified \(x\). They do not, by themselves, give a confidence band with simultaneous coverage over an entire range of \(x\). A parametric bootstrap can propagate parameter uncertainty by repeatedly simulating and refitting under the fitted model, then recalculating the curve. The resulting intervals require the model and bootstrap procedure to be adequate.
Check that numerical integration is sufficiently accurate by increasing the number of quadrature points or Monte Carlo draws and confirming that predictions are stable. Inspect model fit and convergence, the plausibility of the random-effects distribution, and whether clusters provide enough information to estimate slope variation. Curves far beyond the observed covariate range rely heavily on model form and should be labelled as extrapolations. If slopes vary substantially, a single line or a few cluster-specific fits may conceal important heterogeneity.
Keep three prediction targets distinct. The population-averaged curve integrates over the random-effects distribution. A cluster-specific conditional curve uses a specified cluster effect, such as an estimated conditional mode. A zero-random-effects prediction sets \(b=0\) and is conditional on that reference value; it is generally not the marginal curve. Likewise, a GLM without random effects does not represent this cluster-level integration. Each target may answer a useful question, but they should not be reported interchangeably.
Check Your Understanding
Answer each question, then compare with the solutions below.
- In one sentence, distinguish a population-averaged curve from a zero-random-effects prediction.
- For a probit model with independent normal random intercept and slope variances \(0.25\) and \(0.09\), what is \(v(x)\)?
- For the log-link example above, why is the marginal mean at \(x=0\) greater than \(\exp(0.2)\)?
- At \(x=1\) in the logistic numerical example, what is the average probability, and why is it not \(\operatorname{expit}(-0.2)\)?
- Name one diagnostic for numerical integration and one limitation of reporting pointwise intervals along a curve.
Solutions. (1) The population curve averages conditional means over the random-effects distribution, whereas the zero-effects prediction evaluates the conditional mean only at \(b=0\). (2) \(v(x)=0.25+0.09x^2\), because the intercept and slope effects are independent. (3) Integration adds half the random-effect variance to the exponent, so \(m(0)=\exp(0.2+0.125)>\exp(0.2)\). (4) It is \(0.4551\); the inverse link of the average predictor is \(0.4502\), and inverse-link transformation does not generally commute with averaging. (5) Increase quadrature points or Monte Carlo draws and check prediction stability; pointwise intervals do not ensure simultaneous coverage over the full covariate range.