From a Conditional Probability to a Population Probability
In a logistic random-intercept model, an event probability is first specified for a cluster with a particular random intercept. A population-averaged probability answers a different question: what is the event probability at specified covariate values after averaging over the random-intercept distribution? The averaging matters because the inverse-logit function is nonlinear. In general, averaging conditional probabilities is not equivalent to applying the inverse logit to an averaged linear predictor, as discussed in the earlier material on inverse links and averaging.
The target here is the event probability for a randomly selected cluster from the population represented by the model, at fixed covariate values. It is not necessarily the probability for a particular observed cluster, whose random intercept may be informed by that cluster’s data.
The coefficients and variance in this definition are model parameters; in an application they are replaced by estimates. For a binary outcome, the conditional mean is the event probability. Averaging that conditional mean over the random-intercept distribution therefore gives the population event probability by the law of total expectation.
Assumptions and the Numerical Integral
This calculation relies on the model being appropriate for the target population. In particular, conditional on the covariates and random intercept, the outcome is Bernoulli with the stated logistic probability; the random intercept has a normal distribution with mean zero and variance \(\sigma_b^2\); and the covariate values are held fixed during integration. For the usual clustered interpretation, clusters are modeled as independent draws from the random-intercept distribution. Repeated outcomes within a cluster share \(b_j\); conditional independence of those outcomes given the random effect and covariates is a standard additional assumption for the joint model.
The integral is one-dimensional but generally has no elementary closed form. A convenient change of variable writes \(b=\sigma_b Z\), where \(Z\sim N(0,1)\):
This form is useful for numerical integration and simulation. Gauss–Hermite quadrature approximates the expectation by evaluating the integrand at a small number of normal-scale nodes, then combining those values with weights. The approximation can be made more accurate by increasing the number of nodes and checking that the resulting probability stabilizes.
Substitute the specified covariate values into the fixed part of the linear predictor, such as \(\eta(x)=\beta_0+\beta_1x\).
For each possible \(b\), evaluate \(\operatorname{expit}\{\eta(x)+b\}\), then average those probabilities under \(N(0,\sigma_b^2)\).
Describe the result as an event probability for the modeled population at the stated covariate values, rather than as a conditional probability for a cluster with a particular \(b\).
Worked Example: A Marginal Probability at a Specified Covariate Value
Suppose a fitted model has \(\hat\beta_0=0\), \(\hat\beta_1=1\), and \(\hat\sigma_b=1/\sqrt{3}\). We want the population-averaged event probability at \(x=1\). For illustration, use three-point Gauss–Hermite quadrature for the standard normal: nodes \(-\sqrt{3},0,\sqrt{3}\) with weights \(1/6,2/3,1/6\). After multiplying each node by \(\hat\sigma_b\), the random-intercept values are \(-1,0,1\).
The fixed part at \(x=1\) is \(\hat\eta(1)=0+1(1)=1\). The three conditional probabilities are therefore \(\operatorname{expit}(0)=0.5000\), \(\operatorname{expit}(1)=0.7311\), and \(\operatorname{expit}(2)=0.8808\). The quadrature calculation is
Thus the estimated event probability is about 71.75% for a randomly selected cluster at \(x=1\), under this fitted model. The value is an approximation to the integral; it is not the conditional probability at \(b=0\), which here is \(\operatorname{expit}(1)=0.7311\).
Why the Integration Must Be Done on the Probability Scale
The random intercept shifts the linear predictor, but the requested average is an average of probabilities. For a specified \(x\), the operation is to transform each value \(\eta(x)+b\) by the inverse logit and then average over \(b\). Substituting the mean random intercept, \(E(b)=0\), before applying the inverse logit gives \(\operatorname{expit}\{\eta(x)\}\), which is the probability at \(b=0\), not generally the population average.
The source of the difference can be seen by writing the two quantities side by side:
They are not generally equal. This distinction is about the probability at a specified covariate value and the model’s random-intercept distribution; it does not by itself determine the size or direction of any difference for every possible parameter combination.
Worked Example: Integrating at a Negative Covariate Value
Use the same fitted model, \(\hat\beta_0=0\), \(\hat\beta_1=1\), and \(\hat\sigma_b=1/\sqrt{3}\), but now set \(x=-1\). The fixed part is \(\hat\eta(-1)=-1\). At the three quadrature random-intercept values \(-1,0,1\), the linear predictors are \(-2,-1,0\), giving probabilities \(0.1192,0.2689,0.5000\).
Applying the quadrature weights gives
The fitted model therefore gives an event probability of about 28.25% at \(x=-1\), averaged over the random-intercept distribution. The conditional probability at \(b=0\) is instead \(\operatorname{expit}(-1)=0.2689\). The two numbers answer different questions.
Interpretation, Estimation, and Diagnostics
When model parameters are estimated, the plug-in marginal probability is obtained by substituting their estimates into the integral. This gives a point estimate; it does not, by itself, describe uncertainty in that estimate. Uncertainty assessment must account for the estimation of the fixed coefficients and the random-intercept variance. A probability reported without its covariate values, target population, and averaging convention is incomplete.
The marginal probability described here averages over the model’s unconditional random-intercept distribution for a randomly selected cluster. It differs from a cluster-specific prediction that uses an estimated random intercept for an observed cluster. If the scientific target is a particular observed cluster, its information and the corresponding conditional prediction are relevant; integrating over the unconditional distribution instead targets a population of clusters.
A practical numerical check is to recompute the integral with more quadrature nodes. Material changes indicate that the initial approximation may be too coarse, especially when the random-intercept variance is large or the linear predictor places much of the integrand in a region where it changes rapidly. Simulation provides another check: draw many \(b\) values from the fitted normal distribution, calculate the conditional probability for each draw, and average them. The simulation estimate has Monte Carlo error, so it should be compared with an appropriate simulation standard error or repeated with more draws.
Numerical stability cannot repair a poor model. Inspect whether the normal random-intercept distribution and the conditional logistic specification are plausible for the application, and whether the requested covariate values are supported by the observed data. Sparse data, influential clusters, or extrapolation can make a fitted probability unreliable even when the integral is evaluated accurately. Also distinguish uncertainty in the fitted parameters from numerical integration error: increasing quadrature nodes addresses the latter, not the former.
Worked Example: Comparing a Population Probability With the Probability at \(b=0\)
Suppose \(\hat\beta_0=-1\), \(\hat\beta_1=0.5\), and \(\hat\sigma_b=1/\sqrt{3}\), and consider \(x=1\). The fixed part is \(-1+0.5(1)=-0.5\). Using the same three-point quadrature nodes gives random-intercept values \(-1,0,1\), so the linear predictors are \(-1.5,-0.5,0.5\). The corresponding probabilities are \(0.1824,0.3775,0.6225\).
The integrated probability is approximated by
At \(b=0\), the conditional probability is \(\operatorname{expit}(-0.5)=0.3775\). The population-averaged estimate is about 38.58%, while the probability at a zero random intercept is about 37.75%. The first averages across the modeled cluster distribution; the second describes a cluster at the reference value \(b=0\).
Alternatives and Limitations
Quadrature is a deterministic way to evaluate this one-dimensional integral. Monte Carlo integration instead draws random intercepts from the fitted normal distribution and averages the resulting conditional probabilities. It is flexible, but its estimate varies across simulation runs and may require many draws for high precision. Using \(\operatorname{expit}\{\hat\beta_0+\hat\beta_1x\}\) is simpler but answers the \(b=0\) conditional question, not the population-averaged one. A population-averaged model such as one estimated with generalized estimating equations is another approach to marginal effects, but it is a different modeling strategy; its estimates are not obtained by this GLMM integration.
The integral is specific to the assumed logistic random-intercept model and its normal distribution. If the random-effect distribution or linear predictor is misspecified, accurate numerical integration still yields a probability under the wrong model. The calculation also holds covariates fixed: it does not average over a covariate distribution. Averaging over both covariates and random effects would define a different population summary.
Check Your Understanding
For each question, use the definition and three-point quadrature rule above. Treat quadrature results as approximations unless symmetry makes them exact.
- For \(\hat\beta_0=0\), \(\hat\beta_1=1\), and \(\hat\sigma_b=1/\sqrt{3}\), calculate the marginal probability at \(x=0\). How does it compare with the conditional probability at \(b=0\)?
- Using those same parameters, state the three conditional probabilities used to approximate \(\bar p(-1)\), and give their quadrature weights.
- For \(\hat\beta_0=-1\), \(\hat\beta_1=0.5\), and \(\hat\sigma_b=1/\sqrt{3}\), what target does \(\bar p(1)\) describe? Why is it not the probability for every cluster?
- If increasing the number of quadrature nodes changes the estimated probability substantially, what should be checked? Does adding nodes account for uncertainty in the fitted coefficients?
- A researcher substitutes \(b=0\) and reports \(\operatorname{expit}(\hat\beta_0+\hat\beta_1x)\) as the population-averaged probability. Explain the error.
Solutions. (1) At \(x=0\), the three conditional probabilities are \(\operatorname{expit}(-1)=0.2689\), \(0.5000\), and \(\operatorname{expit}(1)=0.7311\), with weights \(1/6,2/3,1/6\). The weighted average is \(0.04482+0.33333+0.12184=0.5000\), matching the conditional probability at \(b=0\); symmetry makes the integrated probability exactly \(0.5\). (2) The probabilities are \(\operatorname{expit}(-2)=0.1192\), \(\operatorname{expit}(-1)=0.2689\), and \(\operatorname{expit}(0)=0.5000\), with weights \(1/6,2/3,1/6\), respectively. (3) It describes the event probability at \(x=1\) for a randomly selected cluster from the population represented by the model, averaged over the random-intercept distribution. Individual clusters have different random intercepts, so their conditional probabilities need not equal that average. (4) Check quadrature convergence and whether the node count is adequate; compare with a simulation estimate if useful. More nodes address numerical integration error, not uncertainty in the fitted coefficients. (5) Substituting \(b=0\) evaluates the probability for a cluster with a zero random intercept. The population-averaged probability requires averaging the inverse-logit probabilities over all random intercepts under their modeled normal distribution.