Marginalizing Interaction Terms

Why an interaction must be marginalized on the response scale

A GLMM interaction describes how predictors combine in the model conditional on the random effects. But many questions concern expected outcomes in a population of clusters, with those random effects averaged over their distribution. With a nonlinear link, integrating over random effects and transforming the linear predictor do not generally commute. Thus the conditional interaction coefficient alone does not determine the size, or even necessarily the direction, of an interaction contrast in marginal response predictions.

The task here is to calculate the response-scale predictions at the predictor combinations of interest, integrate over the random-effects distribution, and form the contrast from those predictions. This is distinct from interpreting the coefficient on the link scale. It also requires specifying the population over which predictions are averaged: for example, a new cluster drawn from the fitted random-effects distribution, and a stated distribution of covariates.

Definition: Let $m(x,a)$ be the marginal mean response at dose $x$ and group $a$, after averaging over the random effects and, if required, a specified target distribution of covariates. For two dose values $x_0,x_1$ and groups $a=0,1$, the response-scale interaction contrast is $$ \Delta_{\mathrm{resp}}= \{m(x_1,1)-m(x_0,1)\} -\{m(x_1,0)-m(x_0,0)\}. $$ It is the difference between the two groups' marginal response changes over the specified dose range. For a binary outcome it is measured in probability units; for a count outcome it is measured in expected counts. It is not, in general, the conditional interaction coefficient.

Integrating a conditional interaction model

Consider a hypothetical binary-outcome GLMM with a random intercept:

Y_{ij}\mid b_j \sim \operatorname{Bernoulli}\{p_{ij}(b_j)\},\qquad \operatorname{logit}\{p_{ij}(b_j)\} =\beta_0+\beta_x x_{ij}+\beta_A A_{ij} +\beta_{xA}x_{ij}A_{ij}+b_j, \qquad b_j\sim N(0,\sigma_b^2).

Here $i$ indexes observations, $j$ indexes clusters, $A$ is a binary group indicator, and $b_j$ is independent of the predictors with the stated normal distribution. The conditional mean is $p_{ij}(b_j)=\operatorname{logit}^{-1}(\beta_0+\beta_xx_{ij}+\beta_AA_{ij}+\beta_{xA}x_{ij}A_{ij}+b_j)$. For a new cluster drawn from this random-effects distribution, its marginal response probability at fixed $x,A$ is

m(x,A)=\int \operatorname{logit}^{-1}\{\beta_0+\beta_xx+\beta_AA+\beta_{xA}xA+b\} \,f_b(b)\,db.

This integral averages probabilities, not linear predictors. In general, $m(x,A)\ne\operatorname{logit}^{-1}\{\beta_0+\beta_xx+\beta_AA+\beta_{xA}xA+E(b)\}$. Even though $E(b)=0$, the latter expression is just the conditional probability at $b=0$, not the marginal probability.

The conditional coefficient $\beta_{xA}$ is a difference in conditional log-odds slopes: the conditional log-odds slope for $A=1$ is $\beta_x+\beta_{xA}$, and for $A=0$ it is $\beta_x$. That statement does not imply that $\beta_{xA}$ equals a difference in marginal probability changes. To obtain that probability contrast, each of the four marginal probabilities must first be calculated.

Numerical integration: A standard-normal three-point Gauss–Hermite approximation uses nodes $0,\pm\sqrt{3}$ and weights $2/3,1/6,1/6$. If $b=\sigma_b Z$, then for a function $g$, $$ E\{g(b)\}\approx \frac{2}{3}g(0)+\frac{1}{6}g(-\sqrt{3}\sigma_b) +\frac{1}{6}g(\sqrt{3}\sigma_b). $$ This is an approximation to the integral, not an identity; more quadrature nodes or simulation may be used when greater numerical accuracy is needed.

Worked marginal probability interaction

Worked Example: A Logistic Random-Intercept Model

Suppose hypothetical fitted parameters are $\beta_0=-1$, $\beta_x=0.5$, $\beta_A=0.4$, $\beta_{xA}=0.6$, and $\sigma_b=0.5$. Let dose $x$ be measured in dose units and compare $x=0$ with $x=2$. Assume the random intercept is normally distributed as specified and independent of dose and group. Use the three-point quadrature above to estimate each marginal probability.

The four conditional linear predictors excluding $b$ are:

$(x,A)=(0,0)$: $-1$.
$(x,A)=(0,1)$: $-1+0.4=-0.6$.
$(x,A)=(2,0)$: $-1+0.5(2)=0$.
$(x,A)=(2,1)$: $-1+0.5(2)+0.4+0.6(2)=1.6$.

The random-effect nodes are $0$ and $\pm\sqrt{3}(0.5)=\pm0.8660$. For instance, at $(0,0)$,

m(0,0)\approx \frac{2}{3}\operatorname{logit}^{-1}(-1) +\frac{1}{6}\operatorname{logit}^{-1}(-1.8660) +\frac{1}{6}\operatorname{logit}^{-1}(-0.1340) \approx 0.2794.

Applying the same formula to each cell gives $m(0,1)\approx0.3618$, $m(2,0)\approx0.5000$, and $m(2,1)\approx0.8209$. In particular, the last predictor is $1.6$, not $1.0$, because the interaction contributes $0.6(2)=1.2$.

The marginal dose change in group $A=0$ is $0.5000-0.2794=0.2206$, a 0.2206 increase in probability over two dose units. In group $A=1$, it is $0.8209-0.3618=0.4591$. Therefore,

\Delta_{\mathrm{resp}}=0.4591-0.2206=0.2385.

The groups' marginal probability changes differ by about 0.2385 over the two-unit dose range, or $0.2385/2=0.1193$ probability units per dose unit as a finite-difference summary. The conditional interaction coefficient is $0.6$ on the log-odds scale per dose unit. Those values differ because they are on different scales and because the logistic transformation and random-effect integration are nonlinear. Neither number is an arithmetic correction to the other.

Standardizing over a target covariate distribution

Integrating over random effects is not always the only averaging needed. If predictions are intended to describe a target population with varying covariates, calculate predictions for each relevant covariate pattern and average them using the target population's weights. The same weights must be used for all four combinations in a contrast if the goal is to compare those combinations in one common target population. The following example treats the stratum-specific probabilities as hypothetical GLMM predictions that have already been integrated over the random effects.

Worked Example: A Covariate-Standardized Interaction Contrast

Suppose a target population has two covariate strata, $z=0$ and $z=1$, with weights 0.70 and 0.30. Hypothetical response-scale predictions from a GLMM, after random-effects integration, are:

Stratum	$m(0,0\mid z)$	$m(1,0\mid z)$	$m(0,1\mid z)$	$m(1,1\mid z)$
$z=0$	0.20	0.35	0.30	0.60
$z=1$	0.50	0.65	0.55	0.80

Standardize each cell using $m(x,A)=0.70m(x,A\mid z=0)+0.30m(x,A\mid z=1)$. Thus:

$m(0,0)=0.70(0.20)+0.30(0.50)=0.14+0.15=0.29$.
$m(1,0)=0.70(0.35)+0.30(0.65)=0.245+0.195=0.44$.
$m(0,1)=0.70(0.30)+0.30(0.55)=0.21+0.165=0.375$.
$m(1,1)=0.70(0.60)+0.30(0.80)=0.42+0.24=0.66$.

The standardized dose changes are $0.44-0.29=0.15$ for $A=0$, and $0.66-0.375=0.285$ for $A=1$. Their difference is $0.285-0.15=0.135$. Equivalently, the four-cell calculation gives $0.66-0.375-0.44+0.29=0.135$. This estimate pertains to the stated 70%/30% target distribution; changing the target weights can change the marginal predictions and their contrast.

When the marginal interaction has a simpler form

The difference between conditional and marginal interactions depends on the outcome distribution and link. A log-link count model with a random intercept provides a useful contrast. Suppose $Y\mid b$ is Poisson with mean $\exp(\eta+b)$, where $b\sim N(0,\sigma_b^2)$ is independent of the predictors. Since a normal random variable has $E\{\exp(b)\}=\exp(\sigma_b^2/2)$,

E(Y\mid x,A)=\exp(\eta)\exp(\sigma_b^2/2).

The common multiplier applies to every predictor combination. Consequently, in this specific model the interaction on the marginal log-mean scale is the same as the conditional interaction coefficient. But response-scale differences in expected counts still depend on the other coefficients and the prediction levels. This special property should not be carried over to logistic models, random-slope models, or settings where the random-effects distribution changes with predictors.

Worked Example: A Log-Link Count Model

Suppose a hypothetical Poisson GLMM has $\eta=\log(2)+\log(1.5)x+\log(2)A+\log(1.25)xA+b$, where $x$ and $A$ are each 0 or 1 and $b\sim N(0,0.25)$. Assume the Poisson conditional mean and normal random intercept are correctly specified and independent. The marginal multiplier from integrating $b$ is $\exp(0.25/2)=\exp(0.125)\approx1.1331$.

The four marginal expected counts are $2(1.1331)=2.2663$ at $(0,0)$, $2(1.5)(1.1331)=3.3994$ at $(1,0)$, $2(2)(1.1331)=4.5326$ at $(0,1)$, and $2(1.5)(2)(1.25)(1.1331)=8.4986$ at $(1,1)$.

The additive response-scale interaction contrast is $(8.4986-4.5326)-(3.3994-2.2663)=3.9660-1.1331=2.8329$ expected counts. In contrast, the conditional log-rate interaction coefficient is $\log(1.25)\approx0.2231$. On the marginal log-mean scale, the difference-in-differences is $\log(1.25)$, because the common factor $\exp(0.125)$ cancels. The additive expected-count contrast is not that log-scale coefficient.

Diagnostics, alternatives, and common mistakes

The integration and contrast are only meaningful relative to a stated model and target population. Check that the link, random-effects distribution, and random-effects dependence assumptions match the fitted GLMM; that the predictor combinations lie in a region where predictions are supported; and that the covariate weights represent the population named in the estimand. For numerical integration, compare quadrature results at increasing numbers of nodes, or check them against simulation, especially when random-effect variances are large or predictions are near the extremes of the response range. A material change as the numerical method is refined signals inadequate integration accuracy, not evidence of a substantive interaction.

The main alternatives answer different questions. A conditional interaction coefficient describes the interaction on the model's link scale at fixed random effects. A marginal response-scale contrast describes a population-averaged difference in response changes. Predictions at $b=0$ describe a cluster with random effect zero, not an average over new clusters. Empirical standardization averages predictions over observed covariates; model-based integration averages over the specified random-effects distribution. Either approach can be appropriate, but the estimand and population must be made explicit.

Common mistake: “The conditional interaction coefficient is 0.6, so the marginal probability interaction is 0.6.” The correction is that 0.6 is a conditional log-odds interaction per dose unit in the logistic example. The marginal probability contrast must be calculated from integrated probabilities and is measured in probability units.
Common mistake: average the linear predictors first and apply the inverse link once. Unless special structure makes this valid, it does not equal averaging response predictions.
Common mistake: use different covariate weights for the four cells without intending to compare different target populations. Use a common target distribution for a directly comparable standardized contrast.

Marginal interactions also have limitations. They depend on the response scale, the dose values or range used for the contrast, the random-effects distribution, and the target covariate distribution. A single contrast summarizes only the specified comparison; it need not describe the interaction at every dose. Point estimates alone do not communicate uncertainty, which requires an appropriate uncertainty procedure for integrated predictions and their joint contrast.

Key takeaway: To marginalize an interaction, first integrate each conditional prediction over the random effects and average over the intended covariate distribution. Then form the response-scale contrast from those predictions. A conditional interaction coefficient is not, in general, that marginal contrast.

Check Your Understanding

Answer each question, then compare with the short solutions below.

In the logistic example, what is the conditional linear predictor at $x=2,A=1$, excluding the random intercept?
Using its reported probabilities, calculate the group $A=0$ marginal probability change from $x=0$ to $x=2$.
Why is $\operatorname{logit}^{-1}(\eta)$ generally not the marginal probability when the random intercept has mean zero?
If the target covariate stratum weights in the standardization example changed, which quantities might change?
In the Poisson example, why does integration preserve the interaction on the marginal log-mean scale, and does that make the additive expected-count contrast equal to $\log(1.25)$?

Solutions. (1) $-1+0.5(2)+0.4+0.6(2)=1.6$. (2) $0.5000-0.2794=0.2206$ probability units. (3) The inverse-logit is nonlinear, so $E\{\operatorname{logit}^{-1}(\eta+b)\}$ generally differs from $\operatorname{logit}^{-1}\{\eta+E(b)\}$. (4) The standardized cell predictions and their interaction contrast can change because they are weighted averages over a different target population. (5) The shared factor $\exp(\sigma_b^2/2)$ cancels in the log-mean difference-in-differences; the additive expected-count contrast is 2.8329 expected counts, not $\log(1.25)\approx0.2231$.

Stratum	\(m(0,0\mid z)\)	\(m(1,0\mid z)\)	\(m(0,1\mid z)\)	\(m(1,1\mid z)\)
\(z=0\)	0.20	0.35	0.30	0.60
\(z=1\)	0.50	0.65	0.55	0.80