Simulation-Based Uncertainty For Marginal Contrasts

Why Simulate Uncertainty in a Marginal Contrast?

A marginal contrast is usually a nonlinear function of several fitted quantities. It may depend on fixed-effect coefficients, the random-effects covariance matrix, integration over the random-effects distribution, and averaging over a target population of covariate profiles. A small change in any fitted quantity can alter the integrated predictions and their difference. Simulation-based uncertainty assessment handles this chain by recalculating the entire estimand at each draw of the fitted parameters.

The key is to distinguish parameter uncertainty from variation among clusters. Random effects describe cluster-to-cluster variation in the model; a marginal prediction integrates over that distribution. Parameter simulation asks how the integrated prediction would vary across plausible fitted parameter values. It does not replace integration with a single random effect or treat the random effect for a particular observed cluster as the target.

Definition: Let $\theta$ collect all fitted model parameters, including fixed effects and the parameters of the random-effects distribution. For target covariate profiles $i=1,\ldots,N$, let $a$ denote a specified exposure or treatment level. A standardized marginal mean is $$ \mu(a;\theta)=\frac{1}{N}\sum_{i=1}^{N}\int g^{-1}\!\left\{x_i(a)^\mathsf{T}\beta+z_i(a)^\mathsf{T}b\right\} \,dF_\theta(b), $$ where $g^{-1}$ is the inverse link, $x_i(a)$ and $z_i(a)$ are the fixed- and random-effects design vectors, and $F_\theta$ is the fitted random-effects distribution. The marginal contrast comparing $a_1$ with $a_0$ is $\Delta(\theta)=\mu(a_1;\theta)-\mu(a_0;\theta)$. Its simulated uncertainty is assessed by repeatedly drawing $\theta$ and recalculating $\Delta(\theta)$.

The definition fixes the target population: the same covariate profiles and weights are used for both exposure levels and, ordinarily, for every parameter draw. The random-effects integral is also part of the estimand. In a logistic random-intercept model, for example, $b\sim N(0,\tau^2)$, and each prediction averages response probabilities over that normal distribution—not over the linear predictors.

Assumptions and the Simulation Procedure

A simulation-based calculation is only as meaningful as the model and parameter-draw distribution it uses. The GLMM must be adequately specified for the intended estimand: the response distribution and link, fixed- and random-effects structure, and random-effects distribution must be defensible. The target covariate profiles and their weights must represent the population about which the contrast is intended to speak. The fitted parameter estimate must be sufficiently identified, and the procedure used to draw parameters must represent its sampling uncertainty reasonably well.

For a likelihood-based fit, a common approach uses a multivariate normal approximation to the sampling distribution of the estimator on an unconstrained parameter scale. Draw the fixed effects and random-effects covariance parameters jointly, using their estimated covariance. For a covariance matrix, an unconstrained parameterization such as its Cholesky factors can be used, with appropriate transformations back to a valid covariance matrix. Drawing variance components independently, or drawing each coefficient from a separate marginal distribution, generally loses dependence and can produce an incorrect contrast distribution.

Specify the estimand.
Choose the two exposure levels, target profiles and weights, response scale, and contrast direction. Record these choices before simulating.

Draw a complete parameter vector.
For each replicate $m$, draw $\theta^{(m)}$ from a justified approximation to the joint sampling distribution of the fitted parameters, or obtain it by refitting bootstrap samples.

Recompute both integrated predictions.
At $\theta^{(m)}$, integrate over the corresponding random-effects distribution for every target profile and exposure level, then apply the specified population weights.

Form the contrast and summarize its draws.
Calculate $\Delta^{(m)}=\mu(a_1;\theta^{(m)})-\mu(a_0;\theta^{(m)})$. Use the distribution of these values to describe and estimate uncertainty.

Integration can be performed by adaptive or Gaussian quadrature, or by Monte Carlo integration. If Monte Carlo integration is used, holding the integration points fixed across parameter draws—often called using common random numbers—can reduce artificial variation due only to changing simulation noise. Increase integration accuracy until the contrast summaries are stable. Parameter draws must not be confused with integration draws: the former represent uncertainty in fitted parameters, while the latter approximate an integral at a chosen parameter value.

Worked Calculations

Worked Example: A Logistic Random-Intercept Contrast

Suppose a hypothetical fitted model has a binary response, a treatment indicator $A$, and a normally distributed random intercept $b\sim N(0,\tau^2)$. The conditional probability is $$ P(Y=1\mid A,b)=\operatorname{logit}^{-1}(\beta_0+\beta_A A+b). $$ For illustration, take the fitted values $\beta_0=-1.00$, $\beta_A=0.80$, and $\tau=0.50$. Compare treatment $A=1$ with $A=0$ for a target population with no additional covariates.

Use the two-point approximation that places equal weight on the standard-normal values $-1$ and $1$. Thus the random-intercept values used here are $-0.50$ and $0.50$. This is a deliberately simple numerical approximation, not an assertion that the normal distribution has only two possible values.

For $A=0$, the two linear predictors are $-1.00-0.50=-1.50$ and $-1.00+0.50=-0.50$. Therefore, $$ \widehat{\mu}(0)=\frac{\operatorname{logit}^{-1}(-1.50)+\operatorname{logit}^{-1}(-0.50)}{2} =\frac{0.1824+0.3775}{2}=0.27995. $$ For $A=1$, the predictors are $-1.00+0.80-0.50=-0.70$ and $-1.00+0.80+0.50=0.30$, giving $$ \widehat{\mu}(1)=\frac{\operatorname{logit}^{-1}(-0.70)+\operatorname{logit}^{-1}(0.30)}{2} =\frac{0.3318+0.5744}{2}=0.45310. $$ The estimated marginal risk difference is $0.45310-0.27995=0.17315$, or about 17.3 percentage points. This is a response-scale contrast averaged over the random-intercept distribution; it is not the conditional treatment coefficient $0.80$.

For uncertainty simulation, each parameter draw would replace all three fitted values, including $\tau$, and both integrated probabilities would be recalculated. The single calculation above is the contrast at the fitted parameter estimate, not an uncertainty estimate by itself.

Worked Example: Estimating Contrast Variability Across Parameter Draws

Continue with the same hypothetical model and the same two-point integration rule. Suppose three joint parameter draws are $(\beta_0,\beta_A,\tau)=(-0.90,0.70,0.60)$, $(-1.00,0.80,0.50)$, and $(-1.10,0.90,0.40)$. These are illustrative draws; three replicates are far too few for a stable uncertainty assessment.

For the first draw, the $A=0$ probabilities are $\operatorname{logit}^{-1}(-1.50)=0.1824$ and $\operatorname{logit}^{-1}(-0.30)=0.4256$, so $\mu^{(1)}(0)=0.3040$. The $A=1$ probabilities are $\operatorname{logit}^{-1}(-0.80)=0.3100$ and $\operatorname{logit}^{-1}(0.40)=0.5987$, so $\mu^{(1)}(1)=0.45435$. Thus $\Delta^{(1)}=0.45435-0.3040=0.1504$.

The second draw gives the fitted-value calculation in the previous example, so $\Delta^{(2)}=0.1732$. For the third draw, $A=0$ has predictors $-1.50$ and $-0.70$, giving $\mu^{(3)}(0)=(0.1824+0.3318)/2=0.2571$. Under $A=1$, the predictors are $-0.60$ and $0.20$, giving $\mu^{(3)}(1)=(0.3543+0.5498)/2=0.45205$. Hence $\Delta^{(3)}=0.1950$.

The mean of these three contrast draws is $(0.1504+0.1732+0.1950)/3=0.1729$. Their sample standard deviation is approximately $$ \sqrt{\frac{(0.1504-0.1729)^2+(0.1732-0.1729)^2+(0.1950-0.1729)^2}{3-1}} \approx 0.0223. $$ This standard deviation summarizes the spread in this tiny illustrative set; it is not a reliable estimate from only three draws. In an analysis, use many more draws and check that summaries stabilize as the number of draws increases.

Worked Example: Standardizing Each Draw Over Two Covariate Profiles

Suppose a hypothetical target population has two equally weighted profiles, $X=0$ and $X=1$. A logistic random-intercept model has $\beta_0=-2.00$, $\beta_A=0.80$, $\beta_X=1.00$, and $\tau=0.40$, with no treatment-by-$X$ interaction. Use the same two-point integration approximation, so the random-intercept values are $-0.40$ and $0.40$.

For profile $X=0$, the $A=0$ probabilities are $\operatorname{logit}^{-1}(-2.40)=0.0832$ and $\operatorname{logit}^{-1}(-1.60)=0.1680$, averaging to $0.1256$. Under $A=1$, the predictors are $-1.60$ and $-0.80$; the probabilities $0.1680$ and $0.3100$ average to $0.2390$. The profile-specific contrast is $0.2390-0.1256=0.1134$.

For $X=1$, the $A=0$ predictors are $-1.40$ and $-0.60$, giving $(0.1978+0.3543)/2=0.2761$. Under $A=1$, they are $-0.60$ and $0.20$, giving $(0.3543+0.5498)/2=0.4521$. This profile's contrast is $0.4521-0.2761=0.1760$. The equally weighted standardized predictions are $\mu(0)=(0.1256+0.2761)/2=0.2009$ and $\mu(1)=(0.2390+0.4521)/2=0.3456$. Their standardized contrast is $0.3456-0.2009=0.1447$, or about 14.5 percentage points.

In a simulation, this entire profile-by-exposure calculation is repeated for every parameter draw. The profile weights remain one-half each here. If the target population instead uses different weights, the estimand changes; changing the weights between draws would add variation that does not represent fitted-parameter uncertainty.

Why the Simulation Estimates Uncertainty

The justification is the same as for estimating the variability of any function of a random vector. If the fitted estimator is $\hat{\theta}$, and its sampling distribution is approximated by a distribution $Q$, then the contrast is a transformed random quantity $\Delta(\theta)$ for $\theta\sim Q$. Drawing $\theta^{(1)},\ldots,\theta^{(M)}$ from $Q$ and recalculating the contrast produces draws from the induced distribution of that function.

In particular, when the draws are independent and $E_Q[\Delta(\theta)^2]$ is finite, the law of large numbers gives $$ \frac{1}{M}\sum_{m=1}^{M}\Delta^{(m)}\longrightarrow E_Q[\Delta(\theta)] \quad\text{and}\quad \frac{1}{M}\sum_{m=1}^{M}\{\Delta^{(m)}-\overline{\Delta}\}^2 \longrightarrow \operatorname{Var}_Q\{\Delta(\theta)\}, $$ as $M$ increases, where $\overline{\Delta}=M^{-1}\sum_m\Delta^{(m)}$. Thus the empirical spread estimates uncertainty induced by the chosen parameter-draw distribution. It does not correct a misspecified GLMM, an unrepresentative target population, or an inaccurate approximation to the estimator’s sampling distribution.

Diagnostics, Alternatives, and Limitations

Check both statistical and computational stability. Confirm that parameter draws respect constraints such as positive variances and positive-definite covariance matrices. Inspect the draw distributions for extreme skewness, implausible values, or concentration near a boundary. Compare the contrast summaries from independent batches or larger values of $M$. Separately increase quadrature accuracy or the number of integration draws: if results change materially, integration error is contaminating the uncertainty calculation.

A multivariate normal approximation can be poor when a variance component is near zero, the likelihood is asymmetric, parameters are weakly identified, or a binary model has separation or sparse data. A cluster-level parametric bootstrap is one alternative: simulate complete data sets under the fitted GLMM, refit each data set, and recalculate the integrated contrast. It can better reflect nonnormal estimator behavior, but may be computationally demanding and inherits assumptions of the fitted model. A Bayesian analysis can instead recalculate the estimand for each joint posterior draw; the resulting spread represents posterior uncertainty and should not be described as a frequentist sampling distribution.

The delta method approximates the variability of a smooth function using local derivatives and the fitted-parameter covariance. It can be efficient and informative when the approximation is adequate, while repeated recalculation naturally accommodates nonlinear integration and contrasts. Neither approach removes the need to define the target estimand carefully. If the target covariate distribution itself is estimated from a sample and inference is meant to include uncertainty in that population distribution, resample or otherwise model that source of uncertainty as well; holding profiles fixed addresses only parameter uncertainty conditional on the chosen target profiles.

Common Pitfall

Common mistake: Draw one random intercept from the fitted distribution for each parameter replicate, calculate a prediction, and call the resulting spread uncertainty in the marginal contrast. That mixes random-effects variation and Monte Carlo integration noise with uncertainty in the fitted parameters. Correction: for each parameter draw, integrate over its random-effects distribution to obtain the marginal predictions, using sufficiently accurate numerical integration. Draw random effects only as numerical integration points when that is the chosen integration method, and control their simulation error.

Key takeaway: Simulation-based uncertainty for a marginal contrast requires a complete joint parameter draw, recalculation of both standardized predictions, and integration over random effects at every draw. The variability of the resulting contrasts reflects the chosen parameter uncertainty approximation, not model error or random-effect variation alone.

Check Your Understanding

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

In a random-intercept logistic model, why is drawing one random intercept per parameter replicate not equivalent to integrating over the random-effects distribution?
A standardized contrast uses three target profiles with weights $0.2$, $0.3$, and $0.5$. What must stay fixed across parameter draws if the goal is parameter uncertainty conditional on this target population?
Why should fixed effects and random-effects covariance parameters be drawn jointly rather than independently?
A contrast summary changes when the number of parameter draws increases from 500 to 2,000, but not when integration accuracy is increased. What should be investigated first?

Solutions. (1) One random intercept is a single draw from the distribution, not an average over it; its variability adds random-effect and integration noise to the parameter uncertainty. Integrate the response prediction over the random-effects distribution at each parameter draw. (2) Keep the three profiles and their weights $0.2$, $0.3$, and $0.5$ fixed, and use them for both exposure levels. (3) The fitted quantities can be correlated; independent draws discard that dependence and may yield an incorrect distribution for the contrast or invalid covariance matrices. (4) Investigate Monte Carlo error from too few parameter draws, and increase the number of draws or compare independent batches.