Marginalizing A GLMM With Random Slopes

Why a Random Slope Changes the Marginal Mean

In a generalized linear mixed model, a random slope allows the relationship between a covariate and the conditional mean to vary across groups. At a particular covariate value, the random intercept and random slope jointly shift the linear predictor. To obtain a population-averaged mean at that value, integrate the conditional mean over the distribution of both random effects. The key simplification is that their combined contribution at a fixed covariate value is a single normally distributed quantity.

This tutorial concerns that integration at one fixed value of the covariate. It does not treat the resulting marginal means as a full curve or as contrasts between covariate values. The earlier material on integrated means and the law of total expectation supplies the general principle: average the conditional mean over the random-effects distribution.

Definition: Consider a GLMM with one covariate \(x\), a random intercept \(b_{0j}\), and a random slope \(b_{1j}\) for group \(j\). Write the conditional linear predictor as \(\eta(x)+b_{0j}+xb_{1j}\), where \(\eta(x)=\beta_0+\beta_1x\). Suppose \((b_{0j},b_{1j})^\mathsf{T}\sim N(\mathbf 0,\Sigma)\), with \(\Sigma=\begin{pmatrix}\tau_{00}&\tau_{01}\\\tau_{01}&\tau_{11}\end{pmatrix}\), and let \(h=g^{-1}\) be the inverse link. The population-averaged mean at fixed \(x\) is \(m(x)=E_{b}\{h[\eta(x)+b_{0j}+xb_{1j}]\}\).

The expectation here is over the model’s random-effects distribution. It is not an average over the observed groups’ fitted random effects, nor is it generally the conditional mean obtained by setting the random effects to zero. The model assumes the random effects have mean zero and the stated joint normal distribution, and that the conditional response mean is specified by the inverse link. The covariance matrix must be positive semidefinite: its diagonal entries are nonnegative and \(\tau_{00}\tau_{11}-\tau_{01}^2\geq0\).

Reduce the Two Random Effects to One

At a fixed \(x\), define \(U(x)=b_{0j}+xb_{1j}\). A linear combination of jointly normal random variables is normal, so \(U(x)\) is normal with mean zero. Its variance follows by expanding the variance of a sum:

v(x)=\operatorname{Var}\{b_{0j}+xb_{1j}\} =\tau_{00}+2x\tau_{01}+x^2\tau_{11}.

Thus the two-dimensional random-effects integral at this covariate value becomes a one-dimensional normal integral. If \(Z\sim N(0,1)\), then \(U(x)\) has the same distribution as \(\sqrt{v(x)}Z\), and

m(x)=\int_{-\infty}^{\infty} h\{\eta(x)+\sqrt{v(x)}z\}\,\phi(z)\,dz =E_Z\left[h\{\eta(x)+\sqrt{v(x)}Z\}\right],

where \(\phi\) is the standard normal density. The calculation is valid when the random-effects distribution and conditional mean model are correctly specified and the integral exists. For common inverse links it does exist; for example, the log-link mean has a finite integral under a normal random effect.

This reduction is useful even when the integral has no elementary closed form. It identifies exactly which part of the random-effects distribution matters at the chosen \(x\): the variance \(v(x)\) of the combined random contribution. The inverse link determines whether that one-dimensional integral can be simplified further.

Worked Example: Finding the Random-Effect Variance at a Fixed Value

Suppose hypothetical fitted random-effects parameters are \(\tau_{00}=0.36\), \(\tau_{11}=0.16\), and \(\tau_{01}=-0.06\). Find the variance and standard deviation of the combined random contribution at \(x=2\).

Substitute into the variance formula:

$$ v(2)=0.36+2(2)(-0.06)+(2^2)(0.16) =0.36-0.24+0.64=0.76. $$

The standard deviation is \(\sqrt{0.76}=0.8718\). Thus at \(x=2\), the random contribution to the linear predictor has distribution \(U(2)\sim N(0,0.76)\), where the second argument is the variance. The covariance matrix is admissible: \(0.36(0.16)-(-0.06)^2=0.0540>0\). This check matters because an invalid covariance matrix cannot describe the assumed joint normal random effects.

Closed-Form Integrals for Two Links

For a log link, \(h(a)=\exp(a)\). If the conditional response mean is \(\exp\{\eta(x)+U(x)\}\), use the moment-generating function of a normal variable: for \(U\sim N(0,v)\), \(E\{\exp(U)\}=\exp(v/2)\). Therefore

m(x)=\exp\{\eta(x)\}E\{\exp[U(x)]\} =\exp\left\{\eta(x)+\frac{v(x)}{2}\right\}.

The extra \(v(x)/2\) is a consequence of averaging on the response scale. It is not an adjustment to the conditional linear predictor for a particular group. For a Poisson model, this gives the population-averaged expected count per observation window when the conditional mean is a count rate or expected count on that same scale.

Worked Example: A Log-Link Mean With a Random Slope

Use the hypothetical covariance values above and suppose \(\eta(x)=-1+0.25x\). For a Poisson log-link model, find the marginal expected count at \(x=2\), per observation window.

First, the fixed part of the linear predictor is \(\eta(2)=-1+0.25(2)=-0.50\), and the combined random-effect variance is \(v(2)=0.76\). The integrated mean is

m(2)=\exp\{-0.50+0.76/2\} =\exp(-0.12)=0.8869.

The population-averaged expected count is therefore about \(0.887\) per observation window. By contrast, setting both random effects to zero gives \(\exp(-0.50)=0.6065\). That is the conditional mean for a group at the random-effects mean, not the integrated population mean. These calculations assume the conditional Poisson mean and log link are appropriate and the stated normal random-effects model holds.

For a probit link, \(h(a)=\Phi(a)\), where \(\Phi\) is the standard normal cumulative distribution function. A useful identity gives a closed form:

E\{\Phi[\eta+U]\}= \Phi\left(\frac{\eta}{\sqrt{1+v}}\right), \qquad U\sim N(0,v).

One justification uses a latent standard normal variable \(\epsilon\), independent of \(U\). Since \(\Phi(\eta+U)=P(\epsilon\leq\eta+U\mid U)\), taking expectations gives \(P(\epsilon-U\leq\eta)\). The difference \(\epsilon-U\) is normal with mean zero and variance \(1+v\), so this probability is \(\Phi\{\eta/\sqrt{1+v}\}\). For a binary probit model this is the integrated event probability, under the conditional probit specification and the assumed normal random effect.

Worked Example: Integrating a Probit Probability

Suppose a hypothetical binary probit model has \(\eta(2)=0.30\) and \(v(2)=0.76\). Find the population-averaged probability of the event at \(x=2\).

Apply the probit identity:

m(2)=\Phi\left(\frac{0.30}{\sqrt{1+0.76}}\right) =\Phi\left(\frac{0.30}{1.3266}\right) =\Phi(0.2261)\approx0.5895.

The model-implied population-averaged event probability at \(x=2\) is about \(0.590\). The conditional probability at zero random effects is \(\Phi(0.30)\approx0.6179\), which is a different estimand. The calculation requires the probit link and normal random contribution; the closed form should not be transferred to a different inverse link without justification.

Logistic-Normal Integration by Quadrature

For the logistic inverse link, \(h(a)=\operatorname{logit}^{-1}(a)=1/(1+\exp(-a))\), the normal integral generally has no elementary closed form. It can be evaluated by numerical integration or approximated by Gaussian quadrature. Quadrature replaces the integral by a weighted sum of inverse-link values at selected normal-distribution nodes.

For example, a five-point Gauss-Hermite rule expressed for \(Z\sim N(0,1)\) uses nodes \(0,\ \pm1.355626,\ \pm2.856970\) and corresponding weights \(0.533333,\ 0.222076\) for each of the first pair, and \(0.011257\) for each of the second pair. The weights sum to approximately one. At fixed \(x\), evaluate the inverse link at \(\eta(x)+\sqrt{v(x)}z\) for each node \(z\), then take the weighted sum. More nodes or adaptive quadrature can be used when a more accurate approximation is required.

Worked Example: Logistic Marginal Probability by Five-Point Quadrature

Suppose a hypothetical logistic random-intercept-and-slope model has \(\eta(x)=0.40\) and \(v(x)=0.64\) at the covariate value of interest. Then \(\sqrt{v(x)}=0.80\). Apply the five-point rule.

The center node gives \(\operatorname{logit}^{-1}(0.40)\approx0.5987\). The two middle nodes give arguments \(0.40\pm0.80(1.355626)\), or approximately \(1.4845\) and \(-0.6845\), whose inverse-logit values are approximately \(0.8152\) and \(0.3353\). The outer nodes give arguments \(0.40\pm0.80(2.856970)\), or approximately \(2.6856\) and \(-1.8856\), with inverse-logit values approximately \(0.9362\) and \(0.1318\).

\begin{aligned} m(x)&\approx 0.533333(0.5987) +0.222076(0.8152+0.3353)\\ &\quad+0.011257(0.9362+0.1318)\\ &=0.3193+0.2555+0.01202\\ &\approx0.5868. \end{aligned}

The quadrature estimate of the population-averaged event probability is \(0.5868\). The plug-in probability at zero random effects is \(0.5987\), so it is not a substitute for this integral. This is a numerical approximation, not an exact identity; additional quadrature points can check whether the reported digits are stable.

Checking the Calculation and Choosing an Approach

The integration method depends on the inverse link. The log link and probit link above yield closed forms under the stated normality assumptions; the logistic link typically requires numerical integration. Simulation is another option: draw many independent \(Z\sim N(0,1)\), compute \(h\{\eta(x)+\sqrt{v(x)}Z\}\) for each draw, and average. Simulation is flexible but introduces Monte Carlo error, whereas quadrature uses fixed, carefully selected nodes and is often more efficient for a one-dimensional integral.

These approaches answer the same model-based question when implemented accurately: what is the mean after averaging over the assumed random-effects distribution? They differ from a conditional prediction that sets \(b_{0j}=b_{1j}=0\), and from a prediction conditional on estimated random effects for a particular observed group. The appropriate choice depends on the estimand, not merely on computational convenience.

Useful checks include confirming that \(v(x)\geq0\), verifying that the covariance matrix is admissible, and checking that numerical integration is stable to increasing the number of nodes. A negative computed \(v(x)\) usually indicates an invalid covariance estimate or an algebra or coding error. If quadrature estimates change materially as nodes are added, the approximation is not yet trustworthy. Simulation estimates should be checked for Monte Carlo precision, for example by repeating the calculation with more draws.

The method also has limits. The integral is only as appropriate as the conditional mean model and random-effects distribution: substantial departures from normality may change the marginal mean. Integration over the specified random-effects distribution does not itself correct misspecification, omitted covariates, or an incorrect link. In addition, the closed forms rely on particular links; they are not generic formulas for every GLMM.

Common Mistakes

Using only the random intercept variance. At fixed \(x\), the random slope also contributes, so use \(v(x)=\tau_{00}+2x\tau_{01}+x^2\tau_{11}\).
Applying the inverse link after averaging the linear predictor. In general, \(E\{h(\eta+U)\}\neq h\{\eta+E(U)\}=h(\eta)\). Integrate the inverse-link mean, not just the linear predictor.
Calling the zero-random-effects prediction population averaged. It describes the conditional mean at \(b_{0j}=b_{1j}=0\); the population-averaged mean integrates over their distribution.
Reporting quadrature as exact. A weighted sum is an approximation. Check numerical stability and report an appropriate number of digits.

Key takeaway: At a fixed covariate value, a normal random intercept and slope combine into one normal random contribution with variance \(v(x)\). The marginal mean is the inverse-link expectation over that distribution; its form is link-specific.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

If \(\tau_{00}=0.25\), \(\tau_{01}=0.02\), and \(\tau_{11}=0.09\), what is \(v(1)\)?
For a log-link model with \(\eta(x)=0\) and \(v(x)=0.50\), what is the marginal mean?
For a logistic model, why is \(h\{\eta(x)\}\) generally not the marginal mean when random-effect variance is positive?
What should be checked if a calculated random-effect variance is negative?

Solutions. (1) \(v(1)=0.25+2(1)(0.02)+1^2(0.09)=0.38\). (2) \(m(x)=\exp\{0+0.50/2\}=\exp(0.25)=1.284\). (3) The inverse-logit function is nonlinear, so averaging \(\operatorname{logit}^{-1}\{\eta+U\}\) is generally not the same as evaluating it at \(U=0\); the normal integral or a numerical approximation is needed. (4) Check the covariance formula, parameter estimates, and positive-semidefinite covariance constraint; a valid variance cannot be negative.