Empirical Standardization With Observed Covariates

From Covariate Records to an Empirical Population Average

A marginal prediction for a GLMM involves two kinds of averaging. First, for a particular covariate record, average the model’s conditional mean over the random-effects distribution. Then average those integrated predictions over the covariate records that define the target population. In empirical standardization, that target distribution is represented by observed covariate records rather than by a single “typical” covariate value.

This distinction matters when the response mean is nonlinear in the linear predictor. In a logistic GLMM, for example, the response probability is a nonlinear function of both the covariates and the random effect. Predicting at the mean covariate values generally does not reproduce the average of the predictions for the actual records.

Definition: Suppose \(m\) observed covariate records define an empirical target population. For record \(r\), let \(x_r(a)\) and \(z_r(a)\) be the fixed- and random-effects design vectors under condition \(a\), and let \(o_r(a)\) be an offset. In a GLMM with \(b\sim N_q(0,D)\), conditional mean \(\mu_r(a,b)=h^{-1}\{o_r(a)+x_r(a)^\mathsf{T}\beta+z_r(a)^\mathsf{T}b\}\), and inverse link \(h^{-1}\), the empirically standardized prediction is \(\Psi_{\mathrm{emp}}(a)=m^{-1}\sum_{r=1}^m E_b\{\mu_r(a,b)\}\). The expectation integrates over the random-effects distribution while holding record \(r\)'s covariates fixed. In practice, parameters are estimated, so the reported prediction substitutes fitted values for \(\beta\) and \(D\).

The record set and the condition \(a\) must be specified. For example, one might use the observed covariates of a study sample and calculate a predicted probability under each of two treatment assignments for every record. Each prediction is integrated over the assumed random-effects distribution; the resulting \(m\) predictions are averaged with equal weight in this definition. This is a model-based population average over the empirical covariate distribution, not a prediction conditional on each record’s estimated cluster effect.

Why Integration Comes Before Averaging Predictions

For each record, the model first gives a conditional mean at a specified random effect \(b\). The marginal mean for that record is its expectation over \(b\). The empirical standardized prediction then averages those \(m\) marginal means. Writing the operations explicitly clarifies what is being averaged:

\Psi_{\mathrm{emp}}(a) = \frac{1}{m}\sum_{r=1}^{m} \int h^{-1}\!\left\{o_r(a)+x_r(a)^\mathsf{T}\beta+z_r(a)^\mathsf{T}b\right\} \,\phi_q(b;0,D)\,db .

Here \(\phi_q(b;0,D)\) is the density of the \(q\)-dimensional random-effects distribution. The covariates remain fixed at record \(r\)'s values inside the integral. The integral usually has no closed form for a logistic inverse link, so it may be evaluated by numerical quadrature or by averaging predictions over random-effect draws. For some links and random-effects structures, analytic integration is available.

The order expresses the estimand, but the finite sum and integral can be interchanged when the relevant expectations exist: both operations are linear. What generally cannot be done is replace all records by their mean covariate values before applying the nonlinear inverse link. In symbols, an average of \(h^{-1}\) applied to record-specific predictors need not equal \(h^{-1}\) applied to an average predictor.

\frac{1}{m}\sum_{r=1}^m E_b\!\left[h^{-1}(\eta_r+z_r^\mathsf{T}b)\right] \ne E_b\!\left[h^{-1}(\bar{\eta}+\bar{z}^\mathsf{T}b)\right] \quad\text{in general.}

The equality can hold in special cases, such as an identity-link model with a linear predictor and a common mean-zero random intercept. It is not a general property of GLMM standardization.

Worked Examples

Worked Example: Standardized Probabilities in a Logistic Random-Intercept Model

Consider two hypothetical covariate records with \(x_1=0\) and \(x_2=2\). A fitted logistic random-intercept model has linear predictor \(-1+0.8x+0.5a+b\), with \(b\sim N(0,0.16)\). Calculate the empirical standardized probabilities under \(a=0\) and \(a=1\), and their difference. 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 scaling by the random-intercept standard deviation \(0.4\), the \(b\) values are approximately \(-0.693,0,0.693\).

At a fixed predictor \(\eta\), the quadrature approximation to the integrated probability is \(q(\eta)=\{g(\eta-0.693)+4g(\eta)+g(\eta+0.693)\}/6\), where \(g(t)=1/(1+e^{-t})\). Under \(a=0\), the two fixed parts are \(-1+0.8(0)=-1\) and \(-1+0.8(2)=0.6\). The corresponding integrated probabilities are approximately \(q(-1)=0.2758\) and \(q(0.6)=0.6407\). Thus the standardized probability is \((0.2758+0.6407)/2=0.4583\).

Under \(a=1\), the fixed parts increase by \(0.5\), giving \(-0.5\) and \(1.1\). The integrated probabilities are approximately \(q(-0.5)=0.3818\) and \(q(1.1)=0.7431\), so the standardized probability is \((0.3818+0.7431)/2=0.5625\). The standardized risk difference is \(0.5625-0.4583=0.1042\), or about 10.42 percentage points.

This calculation integrates over the stated random-intercept distribution for each observed \(x\), then averages the two results. As a check on the contrast, the treatment term adds the same \(0.5\) to each record’s conditional log odds. The standardized risk difference is not itself a conditional odds contrast; it is a difference between two averaged response probabilities.

Worked Example: Analytic Integration for Poisson Rates

Suppose two hypothetical observed records have fitted Poisson log-rate predictors \(\eta_1=-1\) and \(\eta_2=0.6\), before adding a random intercept. The model is \(\log\{\mu_r(b)\}=\eta_r+b\), with \(b\sim N(0,0.25)\). A second condition adds \(0.4\) to each fixed predictor. Find the empirical standardized rate under each condition and the rate difference.

For a normal random intercept, \(E(e^b)=e^{\sigma_b^2/2}\). Therefore \(E_b\{\mu_r(b)\}=e^{\eta_r+\sigma_b^2/2}\). Here \(\sigma_b^2/2=0.125\). Under the first condition, the record-specific marginal rates are \(e^{-1+0.125}=e^{-0.875}\approx0.4169\) and \(e^{0.6+0.125}=e^{0.725}\approx2.065\) rates per exposure unit. The empirical standardized rate is \((0.4169+2.065)/2\approx1.241\) rates per exposure unit.

Under the second condition, the fixed predictors are \(-0.6\) and \(1.0\). The marginal rates are \(e^{-0.6+0.125}=e^{-0.475}\approx0.6219\) and \(e^{1.0+0.125}=e^{1.125}\approx3.080\). Their average is approximately \((0.6219+3.080)/2=1.851\) rates per exposure unit. The standardized rate difference is \(1.851-1.241=0.610\) rates per exposure unit, allowing for rounding.

The random intercept multiplies every record’s rate by \(e^{0.125}\) after integration. The fixed \(0.4\) condition effect also multiplies each rate by \(e^{0.4}\), so the standardized rate ratio is \(e^{0.4}\approx1.492\). The rate difference, unlike the rate ratio, depends on the distribution of rates across the observed records.

Worked Example: When the Mean Covariate Does Give the Same Answer

For two hypothetical records with \(x_1=0\) and \(x_2=2\), suppose a Gaussian identity-link random-intercept model has conditional mean \(E(Y_r\mid b)=4+1.5x_r+b\), with \(b\sim N(0,1)\). Calculate the empirical standardized mean and compare it with the prediction at the mean covariate.

Because \(E(b)=0\), the integrated predictions are \(4+1.5(0)=4\) and \(4+1.5(2)=7\) outcome units. Their empirical average is \((4+7)/2=5.5\) outcome units. The mean covariate is \(\bar{x}=(0+2)/2=1\); the integrated prediction at that value is \(4+1.5(1)=5.5\) outcome units. They agree because the inverse link is the identity, the fixed predictor is linear in \(x\), and the random intercept has mean zero. That agreement should not be assumed for nonlinear response means.

Assumptions, Numerical Checks, and Failure Modes

Empirical standardization is only as credible as the fitted GLMM and the chosen record set. The integration assumes that the specified random-effects distribution, including its covariance matrix \(D\), is suitable for the population being represented. The conditional mean and link must also be appropriate for the outcome. If random slopes are present, their covariates belong in \(z_r(a)\), and the full joint random-effects distribution must be integrated rather than treating each random effect as unrelated.

The covariate records should represent the population to which the standardized result is intended to refer. When evaluating a condition that changes a covariate or treatment value, the model must support predictions at those values. Predictions far outside the observed covariate support can be numerically calculable but heavily dependent on model extrapolation. If a causal interpretation is intended, additional causal assumptions are required; standardization itself does not establish causation.

For numerical integration, check that increasing the quadrature order leaves the predictions stable to a substantively negligible degree. With simulation, use enough random-effect draws to make Monte Carlo variation small relative to the precision being reported, and check that different random-number seeds give similar results. A calculation that integrates only at \(b=0\) has not integrated over random effects. It is a prediction for a cluster with random effect zero, and generally differs from the marginal prediction.

Fix the target records.
List the observed covariate records and state which condition or conditions are being predicted.

Form a prediction for each record.
Use that record’s covariates, offsets, and condition values in the GLMM linear predictor.

Integrate over random effects.
Evaluate the response-scale conditional mean over the fitted random-effects distribution for each record.

Average the integrated predictions.
Take their arithmetic mean for the empirical standardized prediction, and report its response-scale units.

A main alternative is to predict at a representative covariate profile, such as the mean covariate values. That answers a different question and can produce a different result, especially with nonlinear links or interactions. Another alternative is a marginal model such as a GEE, which specifies a population-averaged mean directly rather than obtaining it by integrating a conditional GLMM. These approaches can be useful, but their modeling assumptions and target quantities should not be treated as automatically interchangeable.

Empirical standardization also has limits. Its value depends on which observed records define the empirical covariate distribution; it is not automatically an estimate for every possible target population. The point estimate alone does not describe uncertainty from estimating the model parameters or random-effects distribution. Finally, numerical integration error, model misspecification, and extrapolation are distinct problems: a highly accurate integral cannot correct a poorly specified model or an unsuitable target population.

Common Mistakes

Using \(b=0\) as if it were integration. The zero random effect gives one conditional prediction. The marginal prediction averages over the random-effects distribution.
Replacing records by their average covariates without justification. Under nonlinear inverse links, the prediction at mean covariates usually differs from the mean of record-specific predictions.
Averaging on the wrong scale. For a population-averaged probability or rate, integrate the response-scale conditional mean, then average those response-scale values.
Reporting the result without naming its target. State which observed records were used and what condition was applied; otherwise the empirical population average is not fully described.
Calling a standardized contrast causal by itself. Standardization is a prediction procedure. A causal interpretation needs additional assumptions beyond the GLMM calculation.

Key takeaway: Empirical standardization evaluates the fitted GLMM for each observed covariate record, integrates each response-scale prediction over the random-effects distribution, and averages the integrated predictions. The covariate records define the empirical target distribution; the random-effects integration defines the marginalization over latent cluster variation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In empirical standardization, what is held fixed inside the random-effects integral?
A binary GLMM gives integrated probabilities \(0.20\), \(0.50\), and \(0.80\) for three observed records. What is their empirical standardized probability?
Why can a prediction at the mean covariate differ from the average of record-specific logistic predictions?
In a Poisson log-link random-intercept model with \(b\sim N(0,0.36)\), by what factor does integration multiply the conditional rate at a fixed covariate record?
Does averaging predictions from a GLMM, by itself, establish that a treatment contrast is causal?

Solutions. (1) The record’s covariates, offset, and specified condition are held fixed; the integral averages over the random-effects distribution. (2) The standardized probability is \((0.20+0.50+0.80)/3=0.50\). (3) The logistic inverse link is nonlinear, so applying it to an average predictor need not equal averaging the probabilities obtained from the separate predictors. (4) The factor is \(E(e^b)=e^{0.36/2}=e^{0.18}\approx1.197\). (5) No. A causal interpretation requires additional causal assumptions; standardization alone is a way to calculate model-based averages.