Marginal Treatment Effects With Random Slopes

Why a Random Treatment Slope Changes the Marginal Effect

In a clustered study, treatment effects may differ from one cluster to another. A random treatment slope represents that variation: one cluster may have a stronger conditional response to treatment than another. The population-averaged treatment contrast is not obtained simply by reading the treatment coefficient from the fitted model. Instead, predict each treatment condition for clusters drawn from the specified random-effects distribution, then average the response-scale difference over that distribution.

This distinction matters especially for nonlinear links. Even when the random effects have mean zero, averaging response probabilities is generally not the same as applying the inverse link to an average linear predictor. With a random treatment slope, the treatment and control predictions also involve different combinations of the random effects, so their marginal distributions need not have the same variance.

Model and Marginal Estimand

Definition: Let $A_{ij}\in\{0,1\}$ denote treatment for observation $i$ in cluster $j$, and let $x_{ij}$ be covariates. A random-intercept, random-treatment-slope GLMM has conditional mean $$ g\{E(Y_{ij}\mid A_{ij},x_{ij},b_j)\} = \beta_0+\beta_A A_{ij}+x_{ij}^{\mathsf T}\beta_x+b_{0j}+b_{1j}A_{ij}, $$ where $g$ is the link function and $$ \begin{pmatrix}b_{0j}\\b_{1j}\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}\tau_{00}&\tau_{01}\\\tau_{01}&\tau_{11}\end{pmatrix} \right]. $$ The population-averaged treatment contrast at covariate value $x$, on the response scale, is $$ \Delta(x)=E_b\left[m(1,x,b)-m(0,x,b)\right], $$ where $m(a,x,b)=E(Y\mid A=a,x,b)$, and the expectation is over the joint random-effects distribution.

The covariance matrix must be positive semidefinite; in the usual nonsingular bivariate-normal case, $\tau_{00}>0$, $\tau_{11}>0$, and $\tau_{00}\tau_{11}-\tau_{01}^2>0$. The covariance $\tau_{01}$ matters because the cluster-specific treatment prediction combines the intercept and slope. For treatment $a$, the random contribution is $b_{0j}+a b_{1j}$, with variance $$ v_a=\tau_{00}+2a\tau_{01}+a^2\tau_{11}. $$ Thus $v_0=\tau_{00}$, while $v_1=\tau_{00}+2\tau_{01}+\tau_{11}$.

This definition averages over the same specified cluster population under both treatment conditions. It assumes the random-effects distribution is appropriate for that population and that the conditional mean model is correctly specified. If the random-effects distribution depends on observed covariates, the integration must use the appropriate covariate-conditional distribution. A causal interpretation additionally requires a defensible treatment-assignment design or causal assumptions; model-based marginalization alone does not establish causality.

Probit Integration: A Derivation and Worked Calculation

For a binary outcome with a probit link, the integration over normal random effects has a closed form. This gives a transparent way to see how the random intercept, random slope, and their covariance enter a marginal treatment effect.

Result: Suppose $P(Y=1\mid A=a,x,b)=\Phi(\eta_a(x)+b_0+a b_1)$, where $\Phi$ is the standard normal cumulative distribution function and the random effects have the bivariate-normal distribution above. Then $$ P(Y=1\mid A=a,x) = \Phi\left(\frac{\eta_a(x)}{\sqrt{1+v_a}}\right), $$ where $\eta_a(x)=\beta_0+\beta_Aa+x^{\mathsf T}\beta_x$. The marginal risk difference is the treated marginal probability minus the untreated marginal probability.

To justify the result, introduce an independent standard normal variable $\epsilon$. The conditional probit probability can be written as $P\{\epsilon\leq \eta_a(x)+b_0+a b_1\mid b\}$. Averaging over $b$ is equivalent to evaluating $P\{\epsilon-b_0-a b_1\leq\eta_a(x)\}$. Since $\epsilon-b_0-a b_1$ is normal with mean zero and variance $1+v_a$, standardizing gives the stated marginal probability. This calculation uses the conditional probit model, normal random effects independent of $\epsilon$, and the specified covariance matrix.

Worked Example: Marginal Risk Difference With a Correlated Random Slope

Consider a hypothetical binary-outcome model with no additional covariates, $\beta_0=-1.00$, and $\beta_A=0.70$. Let $\tau_{00}=0.25$, $\tau_{11}=0.36$, and $\tau_{01}=-0.10$. The covariance matrix is valid because its determinant is $0.25(0.36)-(-0.10)^2=0.0800>0$. Find the marginal risk difference.

For control, $v_0=0.25$, so the marginal probability is $$ P(Y=1\mid A=0)=\Phi\left(\frac{-1.00}{\sqrt{1+0.25}}\right) =\Phi(-0.8944)\approx 0.1855. $$

For treatment, the random-effect variance is $$ v_1=0.25+2(-0.10)+0.36=0.410. $$ The fixed linear predictor is $-1.00+0.70=-0.30$, so $$ P(Y=1\mid A=1)=\Phi\left(\frac{-0.30}{\sqrt{1+0.410}}\right) =\Phi(-0.2526)\approx 0.4002. $$

Therefore the population-averaged risk difference is approximately $0.4002-0.1855=0.2147$, or 21.5 percentage points. The conditional treatment coefficient $0.70$ is a probit-scale shift for a cluster with random effects held fixed; it is not this response-scale risk difference. The negative intercept-slope covariance reduces the treated random-effect variance relative to what it would be if the covariance were zero.

Standardizing Over Covariate Profiles

If the target population contains several covariate profiles, first compute an integrated treatment contrast at each profile, then average those contrasts using the target population weights. These weights define whose average effect is being reported. The random-effects integration and the covariate standardization are distinct operations: the former averages over cluster effects, while the latter averages over observed or specified covariate profiles.

Worked Example: Standardizing Over Two Profiles

Continue the hypothetical probit model above, but include a binary covariate $x$ with coefficient $0.80$. Suppose the target population has 70% with $x=0$ and 30% with $x=1$. The random-effects covariance and treatment coefficient are unchanged. Compute the standardized risk difference.

For $x=0$, the marginal probabilities are the ones already calculated: control $0.1855$, treated $0.4002$, and the profile-specific contrast is approximately $0.2147$.

For $x=1$, the control linear predictor is $-1.00+0.80=-0.20$, and the treated linear predictor is $-0.20+0.70=0.50$. Therefore $$ P(Y=1\mid A=0,x=1)=\Phi\left(\frac{-0.20}{\sqrt{1.25}}\right) =\Phi(-0.1789)\approx 0.4290, $$ $$ P(Y=1\mid A=1,x=1)=\Phi\left(\frac{0.50}{\sqrt{1.41}}\right) =\Phi(0.4211)\approx 0.6631. $$ The profile-specific risk difference is approximately $0.6631-0.4290=0.2341$. Using the unrounded model-based contrasts in the weighted average gives $$ \widehat{\Delta}_{G}=0.70(0.2147)+0.30(0.2341)\approx 0.2205. $$ Thus the standardized risk difference is about 0.221, or 22.1 percentage points, for the stated target population. The weights are part of the estimand: a population with a different distribution of $x$ can have a different average contrast.

When the Marginal Contrast Does Equal the Fixed Treatment Coefficient

Random slopes do not automatically make every marginal contrast differ from the fixed treatment coefficient. With an identity link, averaging a linear mean over mean-zero random effects preserves the fixed-effect mean. This special case is useful as a check, but it should not be generalized to nonlinear links.

Worked Example: Identity-Link Outcome With Heterogeneous Effects

Suppose a hypothetical continuous outcome, measured in points, follows $$ E(Y\mid A=a,b)=50+2.40a+b_0+b_1a, $$ where $E(b_0)=E(b_1)=0$, $\operatorname{Var}(b_1)=1.00$, and the random effects have a valid joint distribution. Find the population-averaged treatment contrast and describe the cluster-level variation.

For treatment, integrating over clusters gives $E(Y\mid A=1)=50+2.40+E(b_0)+E(b_1)=52.40$ points. For control, $E(Y\mid A=0)=50+E(b_0)=50.00$ points. The marginal contrast is $52.40-50.00=2.40$ points, equal to $\beta_A$. Yet a cluster’s conditional treatment contrast is $2.40+b_1$ points, whose standard deviation is $\sqrt{1.00}=1.00$ point. The fixed coefficient is the mean of these cluster-specific contrasts, not a claim that every cluster has a 2.40-point effect.

Estimation, Diagnostics, and Alternatives

In practice, fit the GLMM, obtain estimates of the fixed effects and random-effects covariance matrix, and substitute them into the marginalization formula or evaluate the defining integral numerically. The response-scale predictions must be averaged after applying the inverse link. For a probit model the formula above supplies the integral directly. For other links, integration generally requires numerical quadrature or simulation; those computations should use the fitted joint distribution, including the estimated intercept-slope covariance.

Useful checks include whether the random-slope variance is estimable, whether the fitted covariance matrix is positive semidefinite, and whether the model converged to a stable solution. A variance estimate near zero, a correlation near its boundary, large changes under modest model adjustments, or a singular fit can indicate that the data do not support the proposed random-effects structure. Random-slope estimation also needs information about treatment contrasts across enough clusters. If treatment never varies within clusters, a within-cluster treatment slope may be weakly identified or not identified separately from cluster-level structure, depending on the design.

Inspect whether predicted probabilities are sensible across treatment arms and covariate profiles, and whether the assumed random-effects distribution is a plausible approximation. A small number of clusters can make the estimated distribution, especially its slope variance and covariance, unstable. Also verify that the target-population weights represent the intended population rather than silently substituting the sample’s cluster or individual weighting scheme.

Approach	What it summarizes	Important distinction
Conditional treatment coefficient	Change in the linear predictor for a cluster, holding its random effects fixed	Not generally a response-scale population average
Integrated GLMM contrast	Average response-scale contrast over the specified joint random-effects distribution	Depends on the link, covariance structure, and target covariate distribution
Random-intercept-only model	Allows clusters to differ in baseline response but constrains the treatment slope to be common	Cannot represent cluster-to-cluster treatment-effect variation

The integrated GLMM contrast retains a model for cluster heterogeneity and yields a clearly defined target-population average. Its limitation is reliance on the conditional mean model and random-effects distribution; misspecification can distort the marginal contrast. A random-intercept-only model may be simpler, but can understate treatment-effect heterogeneity when slopes truly vary. Neither the conditional coefficient nor the marginal contrast is universally preferable: the appropriate summary depends on whether the scientific question concerns a cluster-conditional effect or an average over a specified population.

Common Pitfall

Common mistake: “The random effects have mean zero, so the marginal treatment effect is just $\beta_A$.” This is valid for the mean contrast in the identity-link example, but not generally for nonlinear links. In a nonlinear model, the inverse link is applied to each conditional predictor before averaging. In addition, treatment changes the random contribution from $b_0$ to $b_0+b_1$, whose variance includes $2\tau_{01}+\tau_{11}$. The correction is to integrate the response-scale predictions under each treatment condition and then subtract them.

A related error is to set $b_0=b_1=0$ and call the resulting contrast population-averaged. That prediction describes a cluster at the random-effects mean, not the average response across the distribution. It coincides with the population average only under special conditions, not by definition.

Check Your Understanding

Answer each question, then compare with the solutions.

In the probit model, what is the variance of the random contribution under treatment, and where does the intercept-slope covariance enter?
Why does a zero-mean random slope not imply that the conditional treatment coefficient is a marginal risk difference?
For an identity-link model with mean $10+3A+b_0+b_1A$ and mean-zero random effects, what is the population-averaged treatment contrast?
What population does a weighted average over covariate profiles describe?

Solutions. (1) Under treatment the random term is $b_0+b_1$, so its variance is $\tau_{00}+2\tau_{01}+\tau_{11}$; the covariance contributes $2\tau_{01}$. (2) The coefficient is a shift on the probit linear-predictor scale, while the risk difference is obtained from two response probabilities after integrating over random effects. (3) The contrast is $3+E(b_1)=3$ outcome units, because the identity link preserves the mean under integration. (4) It describes the population represented by the chosen profile weights; changing those weights can change the standardized contrast.