Why Integrating the Probability Matters
In a binary GLMM with a logit link, the model gives a probability conditional on the random effects. A probability for a new cluster, averaged over the random-effects distribution, is a different quantity: it requires averaging the conditional probability, not just the linear predictor. This integral is generally not available in closed form. A logistic-normal approximation gives a quick estimate, but it should be understood as an approximation to the integral rather than a replacement for the model.
The aim here is to estimate that integrated probability at specified covariate values and random-effects parameters. This is distinct from estimating uncertainty in the prediction. The approximation is most useful when a fast calculation is needed and its accuracy can be checked against numerical integration.
Model and Marginal Probability
Suppose \(Y\) is a binary response and the conditional model is Bernoulli with \(\operatorname{logit}\{P(Y=1\mid b)\}=\eta+z^\mathsf{T}b\). Here \(\eta=x^\mathsf{T}\beta\) is the fixed-effects part of the linear predictor at the covariate values of interest, \(z\) is the random-effects design vector, and \(b\) is a mean-zero Gaussian random-effects vector with covariance matrix \(D\). The fixed-effects design vector is \(x\), and \(\beta\) is the fixed-effects coefficient vector. Any offset can be included in \(\eta\).
For this calculation, assume the specified covariates and model parameters are held fixed, the random-effects distribution is Gaussian as stated, and the conditional event probability uses the logistic link. The random-effects contribution \(U=z^\mathsf{T}b\) is then normal with mean zero and variance \(v=z^\mathsf{T}Dz\). Thus the marginal event probability for a new cluster at these covariates is
Here \(\phi(u;0,v)\) is the normal density with mean zero and variance \(v\). If there is only a random intercept with variance \(\tau^2\), then \(v=\tau^2\). With a random intercept and slope, \(v\) depends on the covariate values through \(z\) and \(D\). The mean-zero assumption is a convention for the random effects; any nonzero mean would need to be included in the linear predictor.
The Logistic-Normal Approximation
For \(v\geq 0\), the approximation used here is
The variance \(v\) is on the linear-predictor scale. It must be computed for the specified covariates and random-effects covariance matrix; it is not generally the variance of a response or a probability. When \(v=0\), the formula reduces to \(\operatorname{expit}(\eta)\), the conditional probability without random-effects variation.
Why the Formula Has This Form
A useful justification represents the logistic function as a probability. Let \(L\) have a standard logistic distribution, independent of \(U\sim N(0,v)\). A standard logistic random variable has variance \(\pi^2/3\). Then
The random variable \(L-U\) has mean zero and variance \(\pi^2/3+v\), because \(L\) and \(U\) are independent and have mean zero. It is not exactly logistic: it is a convolution of logistic and normal distributions. Approximate it by a logistic distribution with scale \(s\). A logistic distribution with scale \(s\) has variance \(s^2\pi^2/3\). Matching variances gives
For a logistic variable with scale \(s\), the probability of being at most \(\eta\) is \(\operatorname{expit}(\eta/s)\). Substituting the matched scale yields the stated approximation. This argument explains the adjustment, but variance matching does not make the convolution exactly logistic. The result is therefore an approximation, not an identity.
Worked Example: A Moderate Linear Predictor and Random Intercept
Suppose a hypothetical binary GLMM has \(\eta=1\) at a selected covariate setting and a random-intercept variance \(v=1\). The conditional probability at a zero random intercept is \(\operatorname{expit}(1)\), but the target is the probability averaged over new-cluster random intercepts. The model assumptions above make the random contribution normal with mean zero and variance 1.
Approximation. Since \(\pi^2\approx9.8696\), the adjusted linear predictor is \[ \frac{1}{\sqrt{1+3(1)/9.8696}} =\frac{1}{\sqrt{1.3040}} \approx 0.8757. \] Thus the approximate marginal probability is \[ \operatorname{expit}(0.8757)=\frac{1}{1+\exp(-0.8757)} \approx 0.7059. \]
Numerical integration for comparison. Direct numerical quadrature of \[ \int \operatorname{expit}(1+u)\phi(u;0,1)\,du \] gives approximately \(0.6967\). The approximation is about \(0.0092\) higher, or just under one percentage point. The comparison is consistent with using the formula as a convenient estimate, not as an exact marginal probability.
Worked Example: A Negative Linear Predictor with Larger Variance
For a second hypothetical setting, let \(\eta=-1\) and let the random-intercept variance be \(v=4\), so the random-intercept standard deviation is \(2\). The fixed part alone gives \(\operatorname{expit}(-1)\approx0.2689\). Because the random-effects distribution is centered at zero and has the stated variance, the integral is still defined by the same model; its value is not obtained by simply retaining that conditional probability.
Approximation. The scale adjustment is \[ \sqrt{1+\frac{3(4)}{\pi^2}} =\sqrt{2.2159} \approx1.4886. \] Therefore \[ p_{\mathrm{marg}} \approx\operatorname{expit}\!\left(\frac{-1}{1.4886}\right) =\operatorname{expit}(-0.6718) \approx0.3381. \]
Numerical integration for comparison. Direct quadrature of \[ \int \operatorname{expit}(-1+u)\phi(u;0,4)\,du \] gives approximately \(0.3480\). The formula is about \(0.0099\) lower. The marginal value exceeds the probability at \(u=0\) in this case; this is possible because averaging a nonlinear response function does not generally preserve its value at the mean of its input.
Worked Example: A Random Slope Makes the Adjustment Covariate-Dependent
Suppose a hypothetical model has a random intercept and random slope, with \[ D= \begin{pmatrix} 0.25&0\\ 0&0.16 \end{pmatrix}, \qquad z= \begin{pmatrix} 1\\ 2 \end{pmatrix}. \] At the selected covariate value, the random-effects contribution has variance \[ v=z^\mathsf{T}Dz =0.25+4(0.16)=0.89. \] Assume the fixed-effects part at this same covariate value is \(\eta=2\). The random-effects contribution is normal with mean zero and variance 0.89 under the stated Gaussian model.
Approximation. The adjusted linear predictor is \[ \frac{2}{\sqrt{1+3(0.89)/\pi^2}} =\frac{2}{\sqrt{1.2705}} \approx1.774. \] Hence \[ p_{\mathrm{marg}}\approx\operatorname{expit}(1.774)\approx0.855. \] Direct numerical quadrature of the corresponding logistic-normal integral gives approximately \(0.852\). The agreement is close for this setting, but the variance calculation shows why the same adjustment cannot automatically be used at every covariate value when random slopes are present.
Checking the Approximation and Choosing an Alternative
A practical check is to evaluate the defining one-dimensional integral numerically at the same \(\eta\) and \(v\), then compare it with the formula. Gaussian quadrature or reliable adaptive quadrature can do this directly. For a random-effects vector, the integral reduces to the one-dimensional form above because \(z^\mathsf{T}b\) is normal; alternatively, numerical integration can be carried out over the full vector. Simulation from the stated random-effects distribution is another option, though it introduces Monte Carlo error.
The following distinctions are useful:
| Approach | What it provides | Main consideration |
|---|---|---|
| Logistic-normal approximation | A quick closed-form estimate from \(\eta\) and \(v\) | Fast, but its error is not known from the formula alone |
| Numerical quadrature | A numerical evaluation of the defining integral | Usually a strong check for a one-dimensional Gaussian contribution |
| Simulation | An average of simulated conditional probabilities | Flexible, but accuracy depends on the number of draws and Monte Carlo variability |
| Conditional prediction at \(b=0\) | \(\operatorname{expit}(\eta)\), the probability at the mean random effect | Not a marginal probability when \(v>0\) |
There is no universal error direction: the approximation may be above or below the numerical integral, as the worked examples illustrate. Check it especially when \(v\) is large, when predictions are near the extremes of the probability scale, or when small probability differences matter. For random slopes, calculate \(v=z^\mathsf{T}Dz\) at each setting rather than reusing a variance from a different setting. A discrepancy with numerical integration is a reason to report the integral instead of relying on the approximation.
The approximation treats the fitted parameters as fixed and addresses only integration over Gaussian random effects. It does not account for uncertainty in estimated coefficients or variance components. It also depends on the assumed Gaussian random-effects distribution and logistic conditional link. If either modeling assumption is inappropriate, a close approximation to the integral under that model does not resolve the model misspecification.
Common Mistake
A common mistake is to average the random effect first, note that \(E(U)=0\), and conclude that the marginal probability is \(\operatorname{expit}(\eta)\). That reasoning applies the nonlinear expit function after averaging its input. The required quantity instead averages the expit of the input: \[ E\{\operatorname{expit}(\eta+U)\} \ne \operatorname{expit}\{E(\eta+U)\} =\operatorname{expit}(\eta) \] in general. The approximation adjusts \(\eta\) using the variance \(v\); direct quadrature evaluates the expectation itself. They coincide when \(v=0\), but need not coincide when random-effects variation is present.
Check Your Understanding
Answer each question using the assumptions and formulas developed above. Short solutions follow the questions.
- If \(\eta=0\), what marginal probability does the approximation give for any finite \(v\)?
- For a random intercept with variance \(v=0.5\), write the scale adjustment used in the formula.
- Why is \(\operatorname{expit}(\eta)\) generally not the integrated probability when \(v>0\)?
- With random intercepts and slopes, what quantity should be recomputed at each covariate setting?
- Does the approximation always overestimate or always underestimate the numerical integral?
Solutions. (1) It gives \(\operatorname{expit}(0)=0.5\), since dividing zero by any positive scale still gives zero. (2) The adjustment is \(\sqrt{1+3(0.5)/\pi^2}\). (3) Because the expit function is nonlinear, averaging \(\operatorname{expit}(\eta+U)\) is generally different from applying expit to the mean input \(\eta+E(U)=\eta\). (4) Recompute \(v=z^\mathsf{T}Dz\), since the random-effects design vector \(z\) changes with the covariate setting. (5) No. The error direction depends on the parameter setting; numerical integration provides a direct comparison.