Uncertainty For Integrated Predictions By The Delta Method

Why an Integrated Prediction Has Parameter Uncertainty

A prediction for a new cluster in a generalized linear mixed model (GLMM) is often obtained by averaging the conditional mean over the model’s random-effects distribution. This produces a population-level prediction rather than a prediction conditional on one particular cluster effect. The integral depends on fitted fixed effects and on the random-effects distribution, so uncertainty in both can affect the prediction.

The delta method provides a local approximation to that uncertainty. It applies a first-order Taylor expansion to the integrated prediction as a function of estimated parameters. The resulting standard error is useful when the estimator is sufficiently regular and its covariance matrix is available. It does not remove uncertainty due to model misspecification, and it does not automatically account for numerical error in evaluating the integral.

Definition: Let \(\boldsymbol{\psi}=(\boldsymbol{\beta}^{\mathsf T},\boldsymbol{\theta}^{\mathsf T})^{\mathsf T}\) collect the fixed-effect parameters \(\boldsymbol{\beta}\) and random-effects distribution parameters \(\boldsymbol{\theta}\). For a specified covariate profile \(x\), an integrated prediction is \(m(x;\boldsymbol{\psi})=E_b\{h(x^{\mathsf T}\boldsymbol{\beta}+z^{\mathsf T}b)\}\), where \(h\) is the conditional mean function, \(b\) has the model’s random-effects distribution indexed by \(\boldsymbol{\theta}\), and \(z\) describes how the random effects enter the linear predictor. The delta-method variance estimate is \(\widehat{\operatorname{Var}}\{m(x;\hat{\boldsymbol{\psi}})\}=\hat{\boldsymbol{g}}^{\mathsf T}\widehat{\operatorname{Var}}(\hat{\boldsymbol{\psi}})\hat{\boldsymbol{g}}\), where \(\hat{\boldsymbol{g}}\) is the gradient of \(m\) evaluated at \(\hat{\boldsymbol{\psi}}\).

The covariance matrix here is the joint covariance matrix, not just separate standard errors for each parameter. Its off-diagonal entries represent covariation between estimates. In particular, fixed-effect and variance-parameter estimates can be correlated; omitting that covariance can either increase or decrease the estimated prediction variance.

Derivation and Conditions

Suppose the estimator \(\hat{\boldsymbol{\psi}}\) is approximately multivariate normal around the true parameter vector \(\boldsymbol{\psi}\), with covariance estimated by \(\widehat{\operatorname{Var}}(\hat{\boldsymbol{\psi}})\). If \(m\) is differentiable near \(\boldsymbol{\psi}\), a first-order Taylor expansion gives the justification for the variance formula.

Expand the prediction.
Near the fitted value, \(m(\hat{\boldsymbol{\psi}})\) is approximately \(m(\boldsymbol{\psi})+\boldsymbol{g}^{\mathsf T}(\hat{\boldsymbol{\psi}}-\boldsymbol{\psi})\), where \(\boldsymbol{g}=\nabla_{\boldsymbol{\psi}}m(\boldsymbol{\psi})\).

Take the variance.
The constant \(m(\boldsymbol{\psi})\) has zero variance, and the variance of the linear term is \(\boldsymbol{g}^{\mathsf T}\operatorname{Var}(\hat{\boldsymbol{\psi}})\boldsymbol{g}\).

Estimate at the fitted model.
Replace the unknown gradient and covariance matrix by estimates evaluated at \(\hat{\boldsymbol{\psi}}\). The square root of the resulting variance is the delta-method standard error.

For a prediction integrated over random effects, differentiating with respect to a variance parameter changes the distribution being averaged over. Under regularity conditions that allow differentiation under the integral sign, the derivative includes that change in distribution, not merely a change in the conditional mean at a fixed random effect. In practice, one can differentiate a closed-form integral or use accurate numerical derivatives of the same integration routine used to calculate the prediction.

The approximation is most defensible when the estimator is approximately normal, the fitted parameter is not near a boundary, the prediction is smooth in the parameters, and the covariance matrix is a reasonable estimate of sampling variability. Variance parameters near zero can violate the usual interior, approximately normal conditions. In that setting, a symmetric Wald interval may be poor even if the standard error calculation is mechanically correct.

Worked Examples

Worked Example: A Probit Random-Intercept Prediction

Suppose a hypothetical probit random-intercept model has conditional success probability \(\Phi(\eta+b)\), where \(\eta\) is the fixed linear predictor and \(b\sim N(0,v)\). For this model, integrating over \(b\) gives \(m(\eta,v)=\Phi\{\eta/\sqrt{1+v}\}\). Suppose the fitted values are \(\hat{\eta}=0.400\) and \(\hat{v}=0.360\), and the estimated covariance matrix for \((\hat{\eta},\hat{v})\) is as shown below.

\widehat{\operatorname{Var}} \begin{pmatrix}\hat{\eta}\\\hat{v}\end{pmatrix} = \begin{pmatrix} 0.0400 & 0.000120\\ 0.000120 & 0.0100 \end{pmatrix}.

The integrated prediction is

\hat{m}=\Phi\left(\frac{0.400}{\sqrt{1.360}}\right) =\Phi(0.3430)\approx 0.6342.

Write \(\phi\) for the standard normal density. Differentiating the prediction with respect to the fitted linear predictor and the random-intercept variance gives

\frac{\partial m}{\partial \eta} =\frac{\phi\{\eta/\sqrt{1+v}\}}{\sqrt{1+v}}, \qquad \frac{\partial m}{\partial v} =-\frac{\eta\,\phi\{\eta/\sqrt{1+v}\}}{2(1+v)^{3/2}}.

At the fitted values, these derivatives are approximately \(0.32255\) and \(-0.0474335\). Thus the estimated variance is

\begin{aligned} \widehat{\operatorname{Var}}(\hat{m}) &=(0.32255)^2(0.0400)+(-0.0474335)^2(0.0100)\\ &\quad+2(0.32255)(-0.0474335)(0.000120)\\ &\approx 0.0041804, \qquad SE(\hat{m})=\sqrt{0.0041804}\approx 0.06466. \end{aligned}

A 95% Wald interval is \(0.6342\pm1.96(0.06466)\), or approximately \((0.5075,0.7609)\). In this hypothetical example, the interval expresses uncertainty about the model-based integrated success probability for this covariate profile. It is not a 95% probability statement about a fixed parameter after the data have been observed.

Worked Example: Averaging Predictions Across Two Covariate Profiles

Suppose the target population gives equal weight to two covariate profiles, with probit linear predictors \(\eta_1=\beta_0\) and \(\eta_2=\beta_0+\beta_1\). Let both predictions be integrated over the same random-intercept variance \(v\). The target average is

M(\beta_0,\beta_1,v) =\frac{1}{2}\left[ \Phi\left\{\frac{\beta_0}{\sqrt{1+v}}\right\} +\Phi\left\{\frac{\beta_0+\beta_1}{\sqrt{1+v}}\right\} \right].

For hypothetical estimates \(\hat{\beta}_0=0\), \(\hat{\beta}_1=0.500\), and \(\hat{v}=0.250\), the two integrated probabilities are \(\Phi(0)=0.5000\) and \(\Phi(0.4472)\approx0.6726\). Therefore \(\hat{M}=(0.5000+0.6726)/2=0.5863\).

At \(\eta_1=0\) and \(\eta_2=0.500\), the derivatives of the individual predictions with respect to their linear predictors are approximately \(0.356825\) and \(0.322868\) (the second rounds to \(0.3229\)). Since \(\beta_0\) shifts both predictors, while \(\beta_1\) shifts only the second, the gradient components are

\frac{\partial M}{\partial\beta_0} =\frac{0.356825+0.322868}{2}\approx0.339847, \qquad \frac{\partial M}{\partial\beta_1} =\frac{0.322868}{2}\approx0.161434.

The first profile’s derivative with respect to \(v\) is zero because its linear predictor is zero. The second profile’s derivative is approximately \(-0.064574\), so \(\partial M/\partial v\approx-0.032287\). Suppose the estimated covariance matrix for \((\hat{\beta}_0,\hat{\beta}_1,\hat{v})\) is diagonal with variances \(0.0100\), \(0.0400\), and \(0.0100\), respectively. Then

\begin{aligned} \widehat{\operatorname{Var}}(\hat{M}) &=(0.339847)^2(0.0100)+(0.161434)^2(0.0400) +(-0.032287)^2(0.0100)\\ &\approx0.0022078, \qquad SE(\hat{M})\approx0.04699. \end{aligned}

The approximate 95% Wald interval is \(0.586320\pm1.96(0.04699)\), or about \((0.4942,0.6784)\). The calculation averages predictions first, then propagates uncertainty for that average. The gradient must reflect how each parameter affects every profile and its assigned weight.

Worked Example: Using a Standard-Deviation Parameter

A model may report the random-intercept standard deviation \(q\) rather than its variance \(v=q^2\). With \(\hat{q}=0.600\), the variance is \(0.360\). The chain rule gives \(\partial m/\partial q=2q(\partial m/\partial v)\). Using the fitted values in the first example,

\frac{\partial m}{\partial q} =2(0.600)(-0.0474335)\approx-0.0569202.

This derivative is with respect to \(q\), not \(v\). To illustrate propagation on the standard-deviation scale, suppose the fitted covariance matrix for \((\hat{\eta},\hat{q})\) is

\begin{pmatrix} 0.0400 & 0.001766\\ 0.001766 & 0.0100 \end{pmatrix}.

The corresponding calculation is

\begin{aligned} \widehat{\operatorname{Var}}(\hat{m}) &=(0.32255)^2(0.0400)+(-0.0569202)^2(0.0100)\\ &\quad+2(0.32255)(-0.0569202)(0.001766)\\ &\approx0.0041291, \qquad SE(\hat{m})\approx0.06426. \end{aligned}

This example uses a covariance matrix reported on the \(q\) scale. A covariance matrix on the \(v\) scale must not be combined directly with a derivative on the \(q\) scale. If one reparameterizes the same fitted estimator, its covariance must also be transformed; the chain rule and a consistently transformed covariance matrix give the same first-order delta-method variance.

Diagnostics, Alternatives, and Limitations

A delta-method standard error should be checked as a calculation, not treated as an automatic software output. The following checks are especially useful for an integrated prediction:

Check the parameter scale. Confirm that each derivative corresponds to the parameterization used by the covariance matrix: for example, variance, standard deviation, or log variance.
Check the full covariance matrix. Include estimated covariances whenever available, and verify that the matrix is symmetric and positive semidefinite up to numerical tolerance.
Check the gradient. Compare analytic derivatives with small, scale-aware finite differences. Large disagreement can indicate a coding error or inaccurate integration.
Check integration accuracy. Increase quadrature accuracy or otherwise verify numerical stability. A small change in the integral does not by itself imply a small standard error.
Check approximation quality. A prediction near zero or one, a variance estimate near its boundary, or a highly nonlinear response can make a symmetric normal approximation misleading.

The delta method is attractive because it is fast and transparent: it identifies which parameter directions contribute to uncertainty and uses the fitted joint covariance directly. Alternatives include likelihood-based intervals and resampling approaches, which can better reflect skewness or boundary behavior but require additional computation and their own assumptions. They answer the same uncertainty question only if they repeat the entire estimation and prediction procedure consistently.

The calculated standard error concerns uncertainty from estimating the parameters in the specified model. It does not include uncertainty about the choice of link, random-effects distribution, covariate distribution, or target population. Nor does it include future outcome variability for a new individual; an integrated mean prediction and a prediction of a future outcome are different quantities. Numerical integration error is also distinct: it concerns approximation of the integral, not sampling variability of the fitted parameters.

Common Pitfall

A common mistake is to propagate only fixed-effect uncertainty because the prediction is reported on the response scale. That reasoning overlooks that the integrated response depends on the random-effects distribution as well. For example, in the probit random-intercept model, changing \(v\) changes the scaling \(\sqrt{1+v}\), so \(\partial m/\partial v\) is generally nonzero even though \(v\) does not appear in the conditional linear predictor.

The correction is to define the integrated estimand first, differentiate it with respect to every estimated parameter that determines it, and use the full joint covariance matrix. A second common error is to average the conditional linear predictors and then transform once; that generally does not equal averaging the response-scale predictions. The estimand’s definition determines which function must be differentiated.

Key takeaway: For an integrated GLMM prediction, the delta-method variance is the gradient of the integrated prediction multiplied on both sides of the joint covariance matrix for fixed effects and random-effects distribution parameters. Correct derivatives, consistent parameter scales, and a credible local normal approximation are all essential.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why can the variance of an integrated prediction depend on a random-intercept variance parameter even when that parameter is absent from the conditional linear predictor?
For a differentiable prediction \(m(\boldsymbol{\psi})\), state the first-order delta-method variance estimate.
If the prediction gradient has two components \(0.20\) and \(-0.10\), and the parameter covariance matrix is \(\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix}\), calculate the delta-method variance and standard error.
Name two situations in which a symmetric Wald interval for an integrated prediction may be unreliable.

Solutions. (1) Integration averages the conditional response over the random-effects distribution, and changing its variance changes that distribution and therefore generally changes the average. (2) The estimate is \(\hat{\boldsymbol{g}}^{\mathsf T}\widehat{\operatorname{Var}}(\hat{\boldsymbol{\psi}})\hat{\boldsymbol{g}}\), with \(\hat{\boldsymbol{g}}\) evaluated at the fitted parameters. (3) The variance is \((0.20)^2(0.04)+(-0.10)^2(0.09)+2(0.20)(-0.10)(0.01)=0.0021\), so the standard error is \(\sqrt{0.0021}\approx0.0458\). (4) Examples include a variance estimate near its zero boundary, a prediction near zero or one, or strong nonlinearity that makes the first-order approximation poor.