Differentiating A Marginal Mean

Why Differentiate an Integrated Mean?

A GLMM specifies a response mean conditional on its random effects. For a covariate $x$, that conditional mean can change with both $x$ and the random effects $b$. The population-averaged mean instead integrates over the random-effects distribution. Its derivative describes the local rate of change in that integrated response mean as $x$ changes, holding other covariates fixed.

The distinction matters because the inverse link is generally nonlinear. Differentiating the conditional mean at an average or “typical” random effect will not usually give the derivative of the integrated mean. The appropriate calculation differentiates each conditional mean and then averages those derivatives, provided the conditions for exchanging differentiation and integration hold.

Definition: Let $m(x,b)=E(Y\mid x,b)=h\{\eta(x,b)\}$ be the conditional response mean, where $h$ is the inverse link and $\eta(x,b)$ is the linear predictor. If the random-effects distribution $F_b$ does not depend on $x$, the integrated mean is $$ M(x)=E_b\{m(x,b)\}=\int m(x,b)\,dF_b(b). $$ When it exists, $M'(x)$ is the marginal derivative: the local change in the integrated response mean per unit increase in $x$.

Here $x$ is a scalar covariate. A partial derivative interpretation applies when other covariates are held fixed. The response-scale derivative has units of the response mean per unit of $x$; for a binary outcome, it is in probability points per unit of $x$, and for a count outcome, it is in expected counts per unit of $x$.

When Can Differentiation Move Inside the Integral?

Differentiating an integral is not automatically valid just because the integrand has a derivative. The derivative must be sufficiently controlled over the random-effects distribution so that changes in $x$ do not produce an uncontrolled contribution from extreme values of $b$.

Theorem: Fix an interior value $x_0$. Suppose that (1) the distribution $F_b$ is fixed as $x$ varies near $x_0$; (2) for almost every $b$ under $F_b$, $m(x,b)$ is differentiable in $x$ on a neighborhood of $x_0$; and (3) on that neighborhood, $\left|\partial m(x,b)/\partial x\right|\leq G(b)$ for almost every $b$, where $E_b\{G(b)\}<\infty$. Then $M$ is differentiable at $x_0$, and $$ M'(x_0)=\int \frac{\partial m(x_0,b)}{\partial x}\,dF_b(b) =E_b\left\{\frac{\partial m(x_0,b)}{\partial x}\right\}. $$

The domination condition is sufficient rather than necessary; other conditions can also justify exchanging differentiation and integration. This version is useful because it names the checks explicitly: differentiability for almost every random effect, an $x$-invariant distribution, and an integrable bound in a neighborhood of the point.

To see why the result holds, write the difference quotient of the integrated mean as an integral of difference quotients. For nonzero $t$, $$ \frac{M(x_0+t)-M(x_0)}{t} =\int \frac{m(x_0+t,b)-m(x_0,b)}{t}\,dF_b(b). $$ For almost every $b$, the integrand converges to $\partial m(x_0,b)/\partial x$ as $t$ approaches zero. By the mean value theorem and the bound in the theorem, the absolute value of each difference quotient is at most $G(b)$. Since $G$ is integrable, the dominated convergence theorem permits taking the limit inside the integral. The resulting limit is the integral of the conditional derivative, which proves the identity.

Applying the Chain Rule

For a differentiable inverse link and predictor, the conditional derivative follows directly from the chain rule:

\frac{\partial m(x,b)}{\partial x} =h'\{\eta(x,b)\}\frac{\partial \eta(x,b)}{\partial x}.

Thus, under the theorem’s conditions,

M'(x)=E_b\left[h'\{\eta(x,b)\}\frac{\partial \eta(x,b)}{\partial x}\right].

For example, if the predictor includes a random slope for $x$, written $\eta(x,b)=\beta_0+\beta_1x+b_0+b_1x$, then $\partial\eta(x,b)/\partial x=\beta_1+b_1$. The conditional derivative is $h'\{\eta(x,b)\}(\beta_1+b_1)$, and the marginal derivative averages this entire product—not its factors separately.

Worked Example: A Poisson Random-Intercept Model

Suppose a hypothetical Poisson GLMM has conditional mean $m(x,b)=\exp(0.3+0.4x+b)$, with random intercept $b\sim N(0,0.25)$. Find the marginal derivative at $x=0$.

The random-effects distribution is fixed in $x$. The conditional mean is differentiable, and its derivative is $0.4\exp(0.3+0.4x+b)$. On any bounded neighborhood of $x=0$, this derivative is bounded by a constant times $\exp(b)$; a normal random variable has finite $E\{\exp(b)\}$. Thus an integrable bound exists.

For a normal random variable with mean zero and variance $0.25$, $E\{\exp(b)\}=\exp(0.25/2)=\exp(0.125)$. Therefore

M(x)=\exp(0.3+0.4x)E\{\exp(b)\} =\exp(0.3+0.4x+0.125),

and

M'(x)=0.4\exp(0.3+0.4x+0.125).

At $x=0$, the exponent is $0.3+0.4(0)+0.125=0.425$. Hence

M'(0)=0.4\exp(0.425)\approx 0.6118.

At $x=0$, the expected count therefore increases locally by about 0.6118 expected counts per one-unit increase in $x$. This is a local marginal rate, not a claim that every unit increase produces exactly the same change over a large range.

Nonlinear Links Require Averaging the Derivative

With a nonlinear inverse link, $h'\{\eta(x,b)\}$ varies with $b$. Consequently, averaging the conditional derivative generally differs from taking the derivative at $b=0$, or from multiplying separately averaged terms. The following example makes the averaging explicit.

Worked Example: A Logistic Model With Two Random-Effect Values

Consider a hypothetical binary-outcome model with conditional probability $m(x,b)=\operatorname{logit}^{-1}(-0.5+1.2x+b)$. Let $b=-1$ or $b=1$, each with probability $1/2$, independent of $x$. Calculate the marginal derivative at $x=0$.

The distribution has finite support and is fixed in $x$, and the logistic mean is differentiable. The derivative is bounded on this finite support, so differentiation and averaging can be exchanged. Since the derivative of the logistic inverse link at linear predictor $\eta$ is $p(1-p)$, where $p=\operatorname{logit}^{-1}(\eta)$, the conditional derivative is $1.2p(1-p)$.

For $b=-1$, the linear predictor at $x=0$ is $-1.5$, so $p\approx0.1824$. The conditional derivative is $1.2(0.1824)(1-0.1824)\approx0.1790$. For $b=1$, the predictor is $0.5$, so $p\approx0.6225$, and the conditional derivative is $1.2(0.6225)(1-0.6225)\approx0.2820$. Averaging gives

M'(0)=\frac{0.1790+0.2820}{2}\approx0.2305.

The marginal probability increases locally by about 0.2305 per one-unit increase in $x$, or about 23.05 probability points per unit. The derivative at $b=0$ would instead use $p=\operatorname{logit}^{-1}(-0.5)\approx0.3775$, giving $1.2p(1-p)\approx0.2820$. That is not the marginal derivative: it evaluates the response at one random-effect value rather than averaging over the specified distribution.

Random Slopes and a Changing Predictor

When $x$ has a random slope, both the linear predictor and its derivative depend on the random effects. The formula remains the same, but the random slope must remain inside the expectation. For some links and distributions, the integral can also be evaluated directly, which provides a useful check.

Worked Example: A Log-Link Model With Correlated Random Slopes

Suppose a hypothetical count model has $\eta(x,b)=0.2+0.3x+b_0+b_1x$, with $(b_0,b_1)$ jointly normal, mean zero, $\operatorname{Var}(b_0)=0.16$, $\operatorname{Var}(b_1)=0.09$, and $\operatorname{Cov}(b_0,b_1)=-0.03$. Find $M'(1)$ for the Poisson log-link mean.

The covariance matrix is valid: its determinant is $0.16(0.09)-(-0.03)^2=0.0135>0$, and its diagonal entries are positive. The distribution is fixed in $x$. The conditional mean is exponential in a normal linear combination, and its derivative has finite expectation on any bounded neighborhood of $x=1$, so the differentiation formula applies.

For a mean-zero normal variable $Z$, $E\{\exp(Z)\}=\exp\{\operatorname{Var}(Z)/2\}$. Here $Z=b_0+b_1x$, with variance

$$ V(x)=0.16+2x(-0.03)+0.09x^2 =0.16-0.06x+0.09x^2. $$

Therefore $M(x)=\exp\{0.2+0.3x+V(x)/2\}$, and differentiating gives

M'(x)=M(x)\{0.3+V'(x)/2\} =M(x)(0.27+0.09x).

At $x=1$, $V(1)=0.19$, so $M(1)=\exp\{0.2+0.3+0.19/2\}=\exp(0.595)\approx1.8130$. The multiplier is $0.27+0.09(1)=0.36$, giving

M'(1)=1.8130(0.36)\approx0.6527.

The expected count increases locally by about 0.6527 expected counts per unit increase in $x$ at $x=1$. The covariance between the random intercept and slope affects $V(x)$, and therefore affects both the integrated mean and its derivative.

Diagnostics, Alternatives, and Limits

For numerical integration or simulation, the target is the average of the conditional derivatives at the specified $x$. A useful independent check is a centered finite difference of integrated means:

\frac{M(x+\epsilon)-M(x-\epsilon)}{2\epsilon}.

This approximation should approach the analytic or numerically integrated derivative as $\epsilon$ decreases, subject to numerical error. If it does not, check the predictor derivative, the inverse-link derivative, the random-effects distribution and covariance, and whether the same integration method and population are used at all covariate values. Monte Carlo estimates can also vary because of simulation error; using the same random-effect draws at nearby values of $x$ can make the finite-difference comparison more stable.

A common mistake is to compute the conditional derivative at the mean random effect and call it the marginal derivative. For nonlinear $h$, the derivative depends on $b$, so in general $E_b\{h'(\eta)\partial\eta/\partial x\}$ is not $h'\{E_b(\eta)\}E_b\{\partial\eta/\partial x\}$. The logistic example shows the discrepancy numerically. Another tempting alternative is to approximate the marginal derivative by finite differences of integrated predictions; this can be helpful as a numerical check, but it does not remove the need to define the integration distribution or control numerical error.

The derivative is model-based and local. It describes change in the specified integrated mean, not necessarily a causal effect. Its interpretation depends on the link, predictor, covariates held fixed, and random-effects distribution used in the integral. The differentiation identity also requires that the distribution be fixed as $x$ changes; if that condition or a suitable integrability condition fails, the formula cannot be used without further justification.

Key takeaway: Under regularity conditions, differentiate the conditional mean with respect to the covariate and then average that derivative over the random-effects distribution. Keep every $b$-dependent term inside the expectation, and interpret the result on the response scale and in the units of $x$.

Check Your Understanding

Try each question, then check the brief solutions below.

In a random-slope model with $\eta(x,b)=\beta_0+\beta_1x+b_0+b_1x$, what is the conditional derivative of $m(x,b)=h\{\eta(x,b)\}$?
Name the three sufficient conditions in the theorem that justify differentiating the integrated mean.
For $M(x)=\exp(0.3+0.4x+0.125)$, calculate $M'(0)$ and state its units if the outcome is a count and $x$ is measured in hours.
Why does a nonlinear inverse link generally make the derivative at $b=0$ different from the marginal derivative?

Solutions. (1) The chain rule gives $h'\{\eta(x,b)\}(\beta_1+b_1)$. (2) The distribution $F_b$ is fixed in $x$; the conditional mean is differentiable in $x$ for almost every $b$ near the point; and an integrable bound controls the absolute conditional derivative there. (3) $M'(0)=0.4\exp(0.3+0.4(0)+0.125)=0.4\exp(0.425)\approx0.6118$ expected counts per hour. (4) A nonlinear link makes $h'\{\eta(x,b)\}$ vary over random effects, so its value at zero need not equal its average over their distribution.