Derivatives of the Approximate Population Mean

Why Differentiate an Approximate Population Mean?

Population-design calculations require sensitivity to population parameters, not only to an individual’s typical prediction. If a fixed effect changes, it can alter the prediction at zero random effects and also alter how predictions curve as random effects vary. That curvature affects the approximate population mean. Differentiating only the typical prediction therefore omits a contribution whenever the second-order Taylor correction depends on the fixed effect.

This tutorial derives derivatives of the second-order approximate population mean. The focus is the mean and its parameter derivatives; the derivative of the approximate population covariance is a separate calculation. The starting point is a smooth prediction function $f_j(\theta,\eta)$ for observation $j$, where $\theta$ is a vector of fixed effects and $\eta$ is a vector of mean-zero random effects with covariance matrix $\Omega$.

Definition: Let $f_j(\theta,\eta)$ be the model prediction for observation $j$. Its second-order Taylor-approximated population mean is $$ m_j(\theta,\Omega) = f_j(\theta,0) + \frac{1}{2}\operatorname{tr}\!\left\{H_j(\theta)\Omega\right\}, $$ where $H_j(\theta)$ is the Hessian matrix of second derivatives of $f_j$ with respect to $\eta$, evaluated at $\eta=0$. The first term is the typical prediction. The trace term is the second-order correction for random-effect curvature.

The approximation follows by expanding $f_j$ around $\eta=0$. The first-order random-effect term has expectation zero because $E(\eta)=0$. The expected quadratic term is $\tfrac12\operatorname{tr}(H_j\Omega)$, since $E(\eta\eta^\mathsf{T})=\Omega$. Thus the formula is exact for a quadratic prediction in $\eta$, and is an approximation for a general smooth prediction when higher-order terms can be neglected.

Differentiating the Mean

Let $q_a$ be any scalar population parameter, such as a fixed effect or an element of a parameterization of $\Omega$. Differentiating the mean requires differentiating every quantity in its formula. In particular, the Hessian can depend on fixed effects, and the covariance matrix can depend on covariance parameters.

\frac{\partial m_j}{\partial q_a} = \frac{\partial f_j(\theta,0)}{\partial q_a} + \frac{1}{2}\operatorname{tr}\!\left\{ \frac{\partial H_j}{\partial q_a}\Omega + H_j\frac{\partial\Omega}{\partial q_a} \right\}.

If $q_a$ is a fixed effect and the covariance parameters are separately parameterized, then $\partial\Omega/\partial q_a=0$. The fixed-effect derivative still includes both the derivative of the typical prediction and the derivative of the random-effect curvature correction. If $q_a$ is a covariance parameter and $f_j$ does not depend on it, the typical prediction and Hessian terms have zero derivative, leaving the covariance contribution.

For clarity, if there are $r$ random effects, the Hessian entries are $H_{j,kl}=\left.\partial^2f_j/\partial\eta_k\partial\eta_l\right|_{\eta=0}$, and the correction is $\tfrac12\sum_{k=1}^r\sum_{l=1}^r H_{j,kl}\Omega_{kl}$. Differentiating a fixed effect inside $H_{j,kl}$ generally requires a mixed third derivative of $f_j$. This does not mean that the prediction itself needs to be approximated to third order; it means the second-order correction changes as the fixed effect changes.

Assumptions: The random effects have mean zero and finite covariance $\Omega$. The prediction is sufficiently differentiable in random effects and population parameters for the Taylor expansion and the displayed derivatives to exist. The approximation is useful when omitted higher-order Taylor terms are small over the random-effect distribution. Differentiating the approximation does not make it exact when that approximation is poor.

Worked Examples

Worked Example: A Fixed Effect Changes Both the Prediction and Curvature

Suppose there is one random effect, $\eta$, with mean zero and variance $\omega^2$, and the hypothetical prediction is $f(\beta,\eta)=\exp(\beta+\beta\eta)$. At $\eta=0$, the typical prediction is $\exp(\beta)$. Since $\partial^2 f/\partial\eta^2=\beta^2\exp(\beta+\beta\eta)$, the Hessian at zero is $\beta^2\exp(\beta)$.

The approximate population mean is therefore $$ m(\beta,\omega^2) = \exp(\beta)\left(1+\frac{1}{2}\beta^2\omega^2\right). $$

Take $\beta=0.400$ and $\omega^2=0.0900$. The typical prediction is $\exp(0.4)=1.4918$, and the curvature correction is $\tfrac12(1.4918)(0.4^2)(0.09)=0.01074$. Thus $m\approx1.5026$.

Differentiating with respect to $\beta$ gives $$ \frac{\partial m}{\partial\beta} = \exp(\beta)\left(1+\frac{1}{2}\beta^2\omega^2+\beta\omega^2\right). $$

At the stated values, the derivative of the typical prediction is $\exp(0.4)=1.4918$. The derivative of the curvature correction is $\exp(0.4)\{0.0072+0.0360\}=0.06445$. Hence $\partial m/\partial\beta\approx1.556$. The extra $0.06445$ would be lost by differentiating only the typical prediction. The derivative with respect to the variance parameter is $\partial m/\partial\omega^2=\tfrac12\beta^2\exp(\beta)=0.1193$ at these values.

Worked Example: A Covariate-Dependent Typical Concentration

Consider a hypothetical one-compartment intravenous bolus prediction at time $t$: $C(t)=D V^{-1}\exp(-CL\,t/V)$. Let $r=WT/70$, and suppose the typical values are $CL=5r^{0.75}$ L/h and $V=50r^{\beta_V}$ L, where $WT$ is in kilograms. Use $D=100$ mg, $WT=100$ kg, $t=4$ h, and $\beta_V=1$. This example first checks the direct typical-prediction derivative; the next example adds its random-effect curvature correction.

Here $r=100/70=1.428571$, so $CL=5(1.428571)^{0.75}\approx6.534$ L/h and $V=50(1.428571)=71.4286$ L. Define $k=CL/V$. Then $k\approx0.09147$ h$^{-1}$, and $$ C(4)=\frac{100}{71.4286}\exp\{-4(0.09147)\} \approx0.9710\ {\rm mg/L}. $$

Because $V=50r^{\beta_V}$, its logarithmic derivative is $\partial\log V/\partial\beta_V=\log r$, while $CL$ does not depend on $\beta_V$. Differentiating the log prediction gives $$ \frac{\partial\log C}{\partial\beta_V} = (-1+kt)\log r, \qquad \frac{\partial C}{\partial\beta_V} = C(-1+kt)\log r. $$

Substituting the values gives $(-1+0.3659)(0.3567)(0.9710)\approx-0.2196$ mg/L per unit of $\beta_V$. Increasing $\beta_V$ slightly at this weight decreases the typical predicted concentration. The units are concentration per unit change in the dimensionless exponent.

Worked Example: Adding the Curvature Derivative to the PK Result

Extend the hypothetical model just used by giving clearance a mean-zero log random effect $\eta$: $CL(\eta)=CL\exp(\eta)$, with variance $\omega^2=0.0400$. Keep the volume typical and fixed. Write $x=kt$, where $k=CL/V$. At $\eta=0$, the prediction is $C$, and its second derivative with respect to $\eta$ is $$ H=C(x^2-x). $$

Indeed, $C(\eta)=D V^{-1}\exp\{-x\exp(\eta)\}$, so $\partial C/\partial\eta=-x\exp(\eta)C(\eta)$ and $\partial^2C/\partial\eta^2=[x^2\exp(2\eta)-x\exp(\eta)]C(\eta)$. Setting $\eta=0$ yields the stated Hessian.

The approximate population mean is $m=C+\tfrac12\omega^2C(x^2-x)$. At the values above, $x\approx0.3659$, so $x^2-x\approx-0.2320$. Therefore $m\approx0.9710+0.0200(0.9710)(-0.2320)=0.9665$ mg/L.

To differentiate, let $L=\log r=0.3567$. From the previous example, $C_{\beta_V}=C(-1+x)L$, and $x_{\beta_V}=-xL$. Differentiating $H=C(x^2-x)$ by the product and chain rules gives $$ H_{\beta_V} = C_{\beta_V}(x^2-x) + C(2x-1)x_{\beta_V} \approx0.0849\ {\rm mg/L}. $$

Thus $m_{\beta_V}=C_{\beta_V}+\tfrac12\omega^2H_{\beta_V} \approx-0.2196+0.0200(0.0849)=-0.2179$ mg/L per unit of $\beta_V$. The curvature term shifts the derivative from the typical-prediction sensitivity, but is relatively small here. Its sign and size depend on the model and the nominal values; it should be calculated, not assumed to vanish.

Derivative Checks and Common Failure Modes

An analytic derivative is most useful when its components can be checked. For a fixed effect, identify separately the derivative at $\eta=0$ and the derivative of the trace correction. For a covariance parameter, check that its effect on $\Omega$ is included. These checks help locate omissions that a plausible-looking final number can conceal.

Omitting fixed-effect dependence in the Hessian. The mistaken argument is that the Hessian concerns random effects, so it cannot change with a fixed effect. In fact, $H_j(\theta)$ is evaluated at zero random effects but can still depend on $\theta$, as in the exponential example.
Differentiating only the typical prediction. This is correct only if the curvature correction has zero derivative with respect to the parameter in question. Otherwise it omits $\tfrac12\operatorname{tr}\{(\partial H_j/\partial q_a)\Omega\}$.
Dropping covariance dependence. If a parameter changes $\Omega$, its contribution $\tfrac12\operatorname{tr}\{H_j(\partial\Omega/\partial q_a)\}$ is part of the mean derivative, even when the prediction function itself does not depend on that parameter.
Using an unreliable local approximation. A Taylor approximation around zero random effects may be poor when predictions are strongly nonlinear over the random-effect distribution. A derivative computed exactly for the approximation can still be a poor approximation to the true population-mean derivative.

A practical numerical check is a centered finite difference of the approximate mean itself: $$ \frac{\partial m_j}{\partial q_a} \approx \frac{m_j(q_a+h)-m_j(q_a-h)}{2h}. $$

Use several suitably small step sizes $h$, recalculating the typical prediction, Hessian, and covariance at each perturbed parameter value. Agreement across step sizes and with the analytic derivative supports the implementation. Large disagreement can indicate a missing chain-rule term, an indexing error, an inappropriate step size, or numerical instability. This check validates differentiation of the stated approximation; it does not validate the approximation against the full random-effects model.

Choosing and Interpreting a Derivative

There are three useful levels of calculation. Differentiating only $f_j(\theta,0)$ is the simplest and may be adequate when the curvature correction is negligible or locally constant. Differentiating the full Taylor approximation retains both the typical and curvature contributions and is the target of this tutorial. Direct integration or simulation of $f_j(\theta,\eta)$, followed by numerical differentiation, can assess the approximation when the random-effect distribution and computation permit it. That approach is more computationally demanding and has its own numerical error; it is a comparison, not an automatic guarantee of accuracy.

Interpret the derivative in the parameter’s scale and the prediction’s units. For example, the PK derivative above is measured in mg/L per unit of $\beta_V$, whereas a derivative with respect to a variance parameter has units of prediction divided by variance units. A derivative is local: it describes the change in the approximate mean for a small parameter change near the specified nominal values. It does not by itself describe a finite change or establish that a design is robust to different nominal values.

Key takeaway: Differentiate the whole approximate mean, not just its typical-prediction term. For fixed effects, account for how they change both the prediction at zero random effects and the random-effect Hessian; for covariance parameters, also account for changes in $\Omega$.

Check Your Understanding

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

For $m=f_0+\tfrac12H\omega^2$, where $f_0$ and $H$ depend on $\beta$, write $\partial m/\partial\beta$ when $\omega^2$ is fixed.
If $H$ does not depend on $\omega^2$, what is $\partial m/\partial\omega^2$?
In the exponential worked example, why does the derivative with respect to $\beta$ exceed the derivative of the typical prediction at the stated positive values?
What does agreement between an analytic derivative and a centered finite difference establish, and what does it not establish?

Solutions. (1) $\partial m/\partial\beta=\partial f_0/\partial\beta+\tfrac12\omega^2\partial H/\partial\beta$; both terms are needed unless the Hessian is locally constant in $\beta$. (2) $\partial m/\partial\omega^2=\tfrac12H$. (3) Here the positive curvature correction increases with $\beta$, adding $\exp(\beta)\{\tfrac12\beta^2\omega^2+\beta\omega^2\}$ to the typical derivative. (4) Agreement supports correct numerical differentiation of the approximate mean at that parameter point. It does not show that the Taylor approximation accurately represents the exact population mean.