Marginal Risk Differences For Binary Outcomes

Why Integrate Before Comparing Risks?

For a binary outcome, the risk at a specified exposure level is the probability that the outcome equals 1. In a generalized linear mixed model, that probability depends on both fixed predictors and random effects. A contrast evaluated at one particular random-effect value describes a conditional comparison. A population-averaged risk difference instead compares event probabilities after averaging over the random-effects distribution. The distinction matters when the scientific question concerns the average risk across clusters rather than the comparison within a cluster having a particular latent effect.

This tutorial considers exposure levels \(A=0\) and \(A=1\), and holds a covariate profile \(x\) fixed. The aim is to integrate the model-predicted event probability at each exposure level over the same random-effects distribution, then subtract the resulting risks. Fixing \(x\) makes the estimand precise without changing the target covariate profile between the two predictions.

For example, \(RD(x)=0.04\) means that the integrated predicted event risk under \(A=1\) is 4 percentage points higher than under \(A=0\), at the specified covariate profile and for the random-effects population represented by \(F_b\). This is a model-based contrast. A causal interpretation additionally requires a causal design or assumptions that justify interpreting the exposure contrast causally; the mixed model alone does not establish those conditions.

Model and Assumptions

A common specification is a logistic random-effects model. Conditional on the random effects and predictors, the outcome is Bernoulli, with

Here \(\operatorname{logit}(p)=\log\{p/(1-p)\}\), \(\beta_0,\beta_A,\beta_x\) are fixed-effect parameters, and \(z(a,x)\) specifies how the random effects enter the linear predictor. The random-effects distribution has mean zero and is commonly multivariate normal, \(b\sim N(0,D)\), where \(D\) is a covariance matrix. The inclusion of \(a\) in \(z(a,x)\) allows, for instance, an exposure random slope as well as a random intercept. The inverse logit is \(\operatorname{expit}(u)=1/(1+e^{-u})\).

The integration requires a specified random-effects distribution and model. The usual interpretation also relies on independent random effects across clusters and on the model's assumptions about how those effects relate to the observed predictors. The conditional Bernoulli model must be appropriate for the binary outcome, and the fixed- and random-effect structure must adequately represent the relevant outcome patterns. Predictions should be made at exposure levels and covariate profiles for which the fitted model is meaningful.

To compare exposure levels as a population-averaged contrast, use the same distribution \(F_b\) for both predictions. This represents averaging each exposure-specific risk over the same modeled population of random effects. If instead different random-effects distributions are used, the subtraction mixes an exposure contrast with a change in the population being averaged over.

Computing the Integrated Difference

For the logistic model, the integrated risks are

and the risk difference is \(RD(x)=\mu(1\mid x)-\mu(0\mid x)\). In general, the logistic function of an averaged linear predictor is not the average of the logistic probabilities. Thus, setting \(b=0\) and computing one probability does not usually produce the integrated risk. The introductory material on marginal probabilities in logistic random-intercept models establishes this distinction; here the focus is the contrast between two such integrated risks.

When the integral has no convenient closed form, it can be approximated numerically. For a discrete approximation with random-effect support points \(b_j\) and weights \(w_j\) summing to 1, compute each predicted probability at each exposure level, take the weighted average separately for each level, and subtract. Numerical quadrature is another common approach for continuous random-effects distributions. The same integration method and random-effects distribution should be used for both exposure levels.

Closed-Form Integration in a Probit Example

Some model and random-effects combinations permit a direct formula. For a probit random-intercept model, suppose \(P(Y=1\mid A=a,x,b)=\Phi(\eta_a+b)\), where \(\Phi\) is the standard normal cumulative distribution function, \(\eta_a\) is the fixed part of the predictor, and \(b\sim N(0,\sigma_b^2)\). The identity \[ E_b\{\Phi(\eta_a+b)\}=\Phi\left(\frac{\eta_a}{\sqrt{1+\sigma_b^2}}\right) \] gives the integrated risk. The assumptions here include the stated probit link and normally distributed random intercept; the formula does not apply automatically to a logistic model.

Inference for an Estimated Risk Difference

In practice, the fixed effects and random-effects distribution are estimated, so the integrated risks and their difference are estimates. Uncertainty must account for estimation of the quantities used in both integrations. A common large-sample approach is the delta method. If \(\theta\) collects the fitted model parameters, including the parameters of \(D\), write the estimated contrast as \(g(\hat{\theta})\). If \(\widehat V\) estimates the covariance matrix of \(\hat{\theta}\), then

The gradient can be evaluated analytically or numerically, with the integrated predictions recomputed as the parameters change. A Wald interval is the estimated difference plus or minus \(1.96\) estimated standard errors for an approximate 95% interval. This approximation relies on the estimator being sufficiently close to a normal sampling distribution and can be unreliable in small samples or when variance parameters are poorly estimated. A parametric bootstrap is one alternative: simulate data under the fitted model, refit it, and recompute the integrated difference for each simulated data set.

Diagnostics, Alternatives, and Limitations

A computed contrast is only as credible as the model and the integration. Check that numerical integration is stable by increasing the number of quadrature points or refining the discrete approximation and confirming that both integrated risks and their difference change negligibly. Also inspect whether the fitted model adequately captures the observed binary-outcome patterns across clusters and predictor values. Strong discrepancies, extreme fitted probabilities, sparse exposure groups, or a random-effects distribution that is poorly supported by the data can make the resulting contrast sensitive to modeling choices.

A frequent error is to calculate both predictions at \(b=0\), subtract them, and call the answer population-averaged. That produces a conditional contrast for a cluster whose random effects equal zero. The correction is to calculate probabilities across the random-effects distribution and average them before interpreting the marginal risks. A second error is to average the linear predictors and then apply the inverse link; in nonlinear models this generally differs from averaging the inverse-link probabilities.

The conditional risk difference at a chosen \(b\) is a legitimate alternative when the scientific question concerns a cluster with that specified random-effect value. A population-averaged estimate from generalized estimating equations is another modeling approach, with a different framework for representing dependence. Neither alternative changes the definition of the integrated GLMM contrast given here. The choice depends on the target question and assumptions, not on which estimate is numerically larger.

Finally, the risk difference is scale-specific and depends on the covariate profile and random-effects distribution used in the integration. It does not by itself describe variation in risks across profiles or clusters, and it does not remove bias from confounding or misspecification. Report the exposure levels, covariate profile, random-effects distribution, integration method, estimate, and uncertainty interval so the estimand can be understood and reproduced.