Worked Workshop: Binary Marginal Effects

From a Fitted Model to Binary Marginal Effects

A logistic random-intercept model describes outcome probabilities conditional on a cluster’s random intercept. A marginal effect for a specified covariate profile instead averages those probabilities over the random-intercept distribution. The calculation has two distinct steps: integrate each profile’s probability, then compare the integrated probabilities on the desired scale.

The target here is a member of a new cluster drawn from the model’s random-intercept distribution, at the stated covariate values. The result is a model-based population average for that profile, not a prediction conditional on a particular observed cluster. Whether a contrast has a causal interpretation requires additional design and causal assumptions; the model calculation alone gives an association.

Definition: Suppose $Y_{ij}$ is a binary response for individual $i$ in cluster $j$, and the fitted random-intercept model is $$ Y_{ij}\mid b_j,\mathbf{x}_{ij}\sim\operatorname{Bernoulli}(p_{ij}),\qquad \operatorname{logit}(p_{ij})=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_j,\qquad b_j\sim N(0,\sigma_b^2). $$ For a specified covariate profile $\mathbf{x}$, its population-integrated probability is $$ \bar p(\mathbf{x})=\int \operatorname{expit}(\mathbf{x}^{\mathsf T}\boldsymbol{\beta}+b)\, \phi(b;0,\sigma_b^2)\,db, $$ where $\operatorname{expit}(u)=1/(1+\exp(-u))$, and $\phi(\cdot;0,\sigma_b^2)$ is the normal density with mean 0 and variance $\sigma_b^2$.

For a binary exposure $A$, hold the other covariates at a stated profile $\mathbf{z}$. The risk difference is $\bar p(1,\mathbf{z})-\bar p(0,\mathbf{z})$. The marginal odds ratio is the odds at $A=1$ divided by the odds at $A=0$, with both odds formed from integrated probabilities. These definitions make the order of operations explicit: average probabilities first, then form odds.

Numerically Integrating the Probability

The logistic-normal integral generally has no simple closed form. One practical method is Gaussian quadrature: evaluate the inverse-logit at selected random-intercept values and take a weighted sum. For a transparent hand calculation, three-point Gauss-Hermite quadrature for $b\sim N(0,\sigma_b^2)$ gives

\bar p(\mathbf{x})\approx \frac{1}{6}\operatorname{expit}(\eta-\sqrt{3}\sigma_b) +\frac{2}{3}\operatorname{expit}(\eta) +\frac{1}{6}\operatorname{expit}(\eta+\sqrt{3}\sigma_b), \qquad \eta=\mathbf{x}^{\mathsf T}\boldsymbol{\beta}.

This is a numerical approximation to the integral, not an identity for the logistic-normal mean. More quadrature points can be used in computation. A result should be stable when the number of points is increased; if it is not, the approximation is not yet adequate.

Worked Example: Risk Difference and Marginal Odds Ratio

Consider a hypothetical fitted model with $\operatorname{logit}(p)=-1.0+0.8A+b$, where $A=0$ denotes unexposed and $A=1$ exposed, and the fitted random-intercept standard deviation is $\sigma_b=0.6$. Calculate the integrated probabilities, risk difference, and marginal odds ratio for this profile.

The model assumptions for this calculation are that the response is Bernoulli conditional on $b$, the random intercept follows the stated normal distribution, and the same random-intercept distribution applies to both exposure profiles. The quadrature nodes are $b=0$ and $b=\pm\sqrt{3}(0.6)\approx\pm1.0392$. For $A=0$, $\eta=-1.0$; the three conditional probabilities at the nodes are approximately 0.1152, 0.2689, and 0.5098.

Thus the integrated probability for $A=0$ is

\bar p(0)\approx \frac{1}{6}(0.115145)+\frac{2}{3}(0.268941)+\frac{1}{6}(0.509806) \approx 0.283453.

For $A=1$, $\eta=-0.2$. The conditional probabilities at the same three nodes are approximately 0.2246, 0.4502, and 0.6983. Therefore,

\bar p(1)\approx \frac{1}{6}(0.224570)+\frac{2}{3}(0.450166)+\frac{1}{6}(0.698303) \approx 0.453923.

The risk difference is $0.453923-0.283453=0.170470$, or about 0.1705. For this profile, the model estimates a 17.05 percentage-point higher population-integrated probability for $A=1$ than for $A=0$.

The marginal odds ratio is

\frac{0.453923/(1-0.453923)}{0.283453/(1-0.283453)} \approx \frac{0.8312}{0.3956} \approx 2.101.

Thus the integrated odds are about 2.101 times as large for $A=1$ as for $A=0$. The contrast is based on probabilities integrated over the random-intercept distribution, not on one cluster’s intercept.

Why the Marginal Odds Ratio Is a Separate Calculation

With no exposure interaction, the conditional odds ratio in the model above is $\exp(0.8)$, about 2.226. It compares the odds at $A=1$ and $A=0$ for a fixed value of $b$. The marginal odds ratio from the worked calculation is about 2.101. It is a different quantity because integrating probabilities and then converting them to odds does not preserve the conditional odds ratio.

This distinction follows directly from the definitions. If $q_a(b)=\operatorname{expit}(\eta_a+b)$, then the marginal probability is $E_b[q_a(b)]$, and the marginal odds are $E_b[q_a(b)]/[1-E_b[q_a(b)]]$. In general, those odds are not the average of the conditional odds $q_a(b)/[1-q_a(b)]$. The inverse-logit and the odds transformation are nonlinear, so averaging and transforming do not commute.

Key result: For a specified covariate profile, calculate $\bar p(0,\mathbf{z})$ and $\bar p(1,\mathbf{z})$ separately by integrating over the random-intercept distribution. Then calculate the risk difference by subtraction or the marginal odds ratio from the odds of those two integrated probabilities. Do not substitute a conditional coefficient for either marginal contrast.

Worked Example: A Marginal Contrast at an Interaction Profile

Suppose a hypothetical fitted model is $\operatorname{logit}(p)=-1.4+0.8A+0.6z-0.5Az+b$, with $\sigma_b=0.7$. Calculate the contrasts at $z=1$. The binary covariate $z$ is held fixed at 1 for both exposure levels; the interaction means the exposure coefficient on the conditional log-odds scale depends on $z$.

At $z=1$, the conditional linear predictor is $-0.8+b$ for $A=0$ and $-0.5+b$ for $A=1$. The quadrature nodes are $0$ and $\pm\sqrt{3}(0.7)\approx\pm1.2124$. The corresponding probabilities for $A=0$, at the negative, zero, and positive nodes, are approximately 0.1179, 0.3100, and 0.6016. Hence

\bar p(0,z=1)\approx \frac{1}{6}(0.1179)+\frac{2}{3}(0.3100)+\frac{1}{6}(0.6016) \approx 0.3266.

For $A=1$, the probabilities at those nodes are approximately 0.1528, 0.3775, and 0.6709, giving

\bar p(1,z=1)\approx \frac{1}{6}(0.1528)+\frac{2}{3}(0.3775)+\frac{1}{6}(0.6709) \approx 0.3890.

The risk difference is approximately $0.3890-0.3266=0.0624$, a 6.24 percentage-point difference. The marginal odds ratio, using the integrated probabilities, is approximately

\frac{0.3890/(1-0.3890)}{0.3266/(1-0.3266)} \approx \frac{0.6367}{0.4847} \approx 1.313.

At $z=1$, the conditional exposure log-odds coefficient is $0.8-0.5=0.3$, so the conditional odds ratio is $\exp(0.3)\approx1.350$. That is not the marginal odds ratio of about 1.313. The difference is not an arithmetic inconsistency: the two odds ratios refer to different comparisons.

A Profile Check: Integration Versus Setting the Intercept to Zero

Setting $b=0$ gives the probability for a cluster whose random intercept is at its distributional mean. It does not generally give the population-integrated probability. The next calculation makes the distinction visible at a different covariate profile.

Worked Example: Comparing a Mean-Intercept Prediction With the Integrated Contrast

Suppose the hypothetical fitted model is $\operatorname{logit}(p)=-1.5+0.4z+0.7A+b$, with $\sigma_b=0.5$. At the specified profile $z=1$, the linear predictors excluding $b$ are $-1.1$ for $A=0$ and $-0.4$ for $A=1$. The three quadrature nodes are $0$ and $\pm0.8660$.

At $A=0$, the conditional probabilities at the negative, zero, and positive nodes are about 0.1227, 0.2497, and 0.4418. Their weighted average is

\bar p(0,z=1)\approx \frac{1}{6}(0.1227)+\frac{2}{3}(0.2497)+\frac{1}{6}(0.4418) \approx 0.2606.

At $A=1$, the three probabilities are approximately 0.2199, 0.4013, and 0.6144, so

\bar p(1,z=1)\approx \frac{1}{6}(0.2199)+\frac{2}{3}(0.4013)+\frac{1}{6}(0.6144) \approx 0.4066.

The population-integrated risk difference is about $0.4066-0.2606=0.1460$. The population-integrated odds ratio is approximately $[0.4066/(1-0.4066)]/[0.2606/(1-0.2606)]\approx1.944$. In contrast, setting $b=0$ gives probabilities $\operatorname{expit}(-1.1)\approx0.2497$ and $\operatorname{expit}(-0.4)\approx0.4013$, whose difference is about 0.1516. These are mean-intercept conditional predictions, not integrated population probabilities.

Assumptions, Diagnostics, and Alternatives

The calculation is only as relevant as its model and target. It relies on the fitted logistic mean structure and on the specified random-intercept distribution being appropriate for the population being averaged over. It also assumes that the covariate profile is meaningful and represented adequately by the fitted model. Conditional Bernoulli outcomes and the stated cluster structure are part of the model; a wrong link, omitted interaction, or unsuitable random-effect distribution can make the computed contrast misleading even when the integration is numerically accurate.

Several checks help separate numerical problems from modeling problems:

Check quadrature convergence. Recalculate with more integration points. Material movement in the probabilities or contrasts signals that the numerical approximation needs refinement.
Check the target distribution. Confirm that the random-effect distribution used for integration represents the intended new-cluster population. Integrating over a different distribution changes the estimand.
Check the profile and coding. Verify exposure coding, interaction terms, and units for continuous covariates. A coefficient or profile entered on the wrong scale changes every linear predictor.
Check model fit and scope. Examine whether the model captures important structure and whether the profile is within a defensible range of the observed covariates. Integration does not repair misspecification or unsupported extrapolation.
Keep estimation uncertainty separate from integration. The quadrature calculation evaluates a fitted model at its parameter estimates. A point estimate alone does not account for uncertainty in those estimates.

Numerical quadrature is usually efficient for a one-dimensional random intercept. Simulation can approximate the same integral by drawing many $b$ values from the fitted normal distribution and averaging the resulting probabilities; this is useful in more complex settings but introduces simulation error. A logistic-normal approximation can be faster, but should be checked against more accurate integration when precision matters. Conditional coefficients answer cluster-specific questions, whereas integrated contrasts answer population-averaged questions. Neither scale is universally preferable; the scientific target determines which to report.

A frequent mistake is to calculate $\exp(\hat\beta_A)$ and label it the marginal odds ratio, or to compute probabilities at $b=0$ and label them population averages. The correction is procedural: for each exposure value, integrate the response probability over the specified random-effect distribution, then form the requested contrast. State the covariate profile and target population so the reported effect has a clear interpretation.

Check Your Understanding

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

In a logistic random-intercept model, what does $\bar p(\mathbf{x})$ average over?
For two exposure profiles, should you average the conditional odds and take their ratio, or first integrate the probabilities and then form the odds?
Why can $\exp(\hat\beta_A)$ differ from a marginal odds ratio even when the model has no exposure interaction?
At a specified profile, the integrated probabilities are 0.20 for $A=0$ and 0.35 for $A=1$. Calculate the risk difference and marginal odds ratio.

Solutions. (1) It averages the model’s conditional response probability over the random-intercept distribution for the stated covariate profile. (2) Integrate each profile’s probabilities first, then form the odds and their ratio; the odds transformation is nonlinear. (3) $\exp(\hat\beta_A)$ compares conditional odds at a fixed random intercept, whereas the marginal odds ratio is computed from integrated probabilities. (4) The risk difference is $0.35-0.20=0.15$, or 15 percentage points. The marginal odds ratio is $[0.35/0.65]/[0.20/0.80]\approx0.5385/0.2500=2.154$. The odds at $A=1$ are therefore about 2.154 times those at $A=0$ for this profile.