Derivatives of the Approximate Population Covariance

Why Differentiate the Covariance?

In a nonlinear mixed-effects model, the approximate population covariance describes how observations vary across subjects and within a subject under a specified sampling design. Its parameters can affect that covariance in several ways: they may change the model sensitivities to random effects, the random-effect covariance itself, or the residual-error covariance. A derivative that includes only one of these routes can therefore omit an important part of how the model predicts variability.

This tutorial differentiates the first-order approximate covariance \(V=G\Omega G^\mathsf{T}+R\). The design and observation times are held fixed. The parameter being differentiated may be a fixed-effect parameter, a random-effect variance or covariance parameter, or a residual-error parameter. If a parameter affects more than one component, all of its effects must be included.

Definition: Let \(m\) be the number of observations per subject and \(q\) the number of random effects. The sensitivity matrix \(G\) is an \(m\times q\) matrix, \(\Omega\) is the \(q\times q\) random-effect covariance matrix, and \(R\) is the \(m\times m\) residual-error covariance matrix. The approximate population covariance is the \(m\times m\) matrix \(V=G\Omega G^\mathsf{T}+R\). For a scalar model parameter \(\psi\), its covariance derivative is the elementwise derivative \(\partial V/\partial\psi\), evaluated with the design held fixed.

This expression follows from a first-order expansion of the individual model around zero random effects, with random effects having mean zero and covariance \(\Omega\), and residual error having covariance \(R\). The derivative calculations below require the components to be differentiable with respect to the parameter at the point of evaluation. The covariance matrices must be valid at that point: in particular, \(\Omega\) and \(R\) are symmetric positive semidefinite. These conditions do not require their derivatives to be positive semidefinite.

The Product Rule for the Covariance

Write \(G_\psi=\partial G/\partial\psi\), \(\Omega_\psi=\partial\Omega/\partial\psi\), and \(R_\psi=\partial R/\partial\psi\). Differentiating each factor in \(G\Omega G^\mathsf{T}\) in turn gives

\frac{\partial V}{\partial\psi} = G_\psi\Omega G^\mathsf{T} + G\Omega_\psi G^\mathsf{T} + G\Omega G_\psi^\mathsf{T} + R_\psi.

The first and third terms arise because the sensitivity matrix occurs on both sides of \(\Omega\). The middle term captures changes in random-effect variability, and the final term captures changes in residual variability. When \(\psi\) affects only one component, the other component derivatives are zero; when it affects several, the terms add.

Key result: At a fixed design, the derivative of the approximate population covariance is the sum of its derivatives through the left sensitivity matrix, the random-effect covariance, the right sensitivity matrix, and the residual covariance. A parameter-dependent component must not be treated as constant merely because it is not the component named for that parameter.

The derivative \(G_\psi\) may require more than the first derivatives used to construct \(G\). If \(G\) is the derivative of the individual prediction \(f\) with respect to random effects \(\eta\), evaluated at \(\eta=0\), then \(G_\psi\) involves mixed derivatives of \(f\) with respect to \(\psi\) and \(\eta\), where they exist. This dependence is one reason covariance derivatives can be more involved than mean derivatives.

Why the Formula Holds

A short differential argument makes the three matrix-product terms explicit. For an infinitesimal change in \(\psi\), the product rule for matrices gives the change in \(G\Omega G^\mathsf{T}\) by changing one factor at a time. Since the transpose operation commutes with differentiation, the change in the last factor is the transpose of the change in \(G\).

\begin{aligned} d(G\Omega G^\mathsf{T}) &=(dG)\Omega G^\mathsf{T} +G(d\Omega)G^\mathsf{T} +G\Omega(dG)^\mathsf{T},\\ dV&=d(G\Omega G^\mathsf{T})+dR. \end{aligned}

For a scalar parameter, \(dG=G_\psi\,d\psi\), \(d\Omega=\Omega_\psi\,d\psi\), and \(dR=R_\psi\,d\psi\). Dividing the resulting expression for \(dV\) by \(d\psi\) yields the stated derivative. This derivation also shows why dropping either sensitivity term is generally incorrect: each corresponds to a distinct occurrence of \(G\) in the covariance product.

If \(\Omega\) and \(R\) are symmetric, then \(\partial V/\partial\psi\) is symmetric. The first and third terms are transposes of one another, and the derivatives of symmetric \(\Omega\) and \(R\) are symmetric wherever the parameterization preserves symmetry. This provides a useful check on an analytic or numerical calculation.

Worked Example: A Parameter Changes Sensitivities

Suppose there are two observations and one random effect, with \(\psi\) dimensionless and \(G(\psi)=(\psi,\,2\psi)^\mathsf{T}\), \(\Omega=0.25\), and \(R=\operatorname{diag}(0.10,0.20)\). Calculate \(V\) and \(\partial V/\partial\psi\) at \(\psi=2\). The random-effect variance is positive and \(R\) has positive diagonal entries, so both covariance components are valid.

Here \(G=(2,4)^\mathsf{T}\), \(G_\psi=(1,2)^\mathsf{T}\), \(\Omega_\psi=0\), and \(R_\psi=0\). First,

G\Omega G^\mathsf{T} =0.25 \begin{pmatrix}2\\4\end{pmatrix} \begin{pmatrix}2&4\end{pmatrix} = \begin{pmatrix}1&2\\2&4\end{pmatrix}.

Thus \(V=G\Omega G^\mathsf{T}+R=\begin{pmatrix}1.10&2.00\\2.00&4.20\end{pmatrix}\), in squared observation units. Its derivative is

\begin{aligned} \frac{\partial V}{\partial\psi} &=G_\psi\Omega G^\mathsf{T}+G\Omega G_\psi^\mathsf{T}\\ &=0.25\left[ \begin{pmatrix}1\\2\end{pmatrix}\begin{pmatrix}2&4\end{pmatrix} + \begin{pmatrix}2\\4\end{pmatrix}\begin{pmatrix}1&2\end{pmatrix} \right]\\ &=0.25\left[ \begin{pmatrix}2&4\\4&8\end{pmatrix} + \begin{pmatrix}2&4\\4&8\end{pmatrix} \right] = \begin{pmatrix}1&2\\2&4\end{pmatrix}. \end{aligned}

The entries of this derivative are changes in covariance per unit change in \(\psi\). The result is symmetric, as required. The residual covariance contributes to \(V\), but not to this derivative because \(R\) is constant in \(\psi\).

Parameters in the Random-Effect Covariance

A parameter can change \(V\) even when the sensitivities remain fixed. For example, a covariance parameter may change the association between two random effects. The following example uses an unconstrained covariance entry \(c\); its value is chosen so that \(\Omega\) is positive definite at the evaluation point.

Worked Example: A Random-Effect Covariance Changes

Let \(G=\begin{pmatrix}1&1\\1&0\end{pmatrix}\), \(\Omega=\begin{pmatrix}0.40&c\\c&0.25\end{pmatrix}\), and \(R=\operatorname{diag}(0.10,0.20)\). Evaluate \(V\) and its derivative with respect to \(c\) at \(c=0.10\). The determinant of \(\Omega\) is \(0.40(0.25)-0.10^2=0.09>0\), and its diagonal entries are positive, so \(\Omega\) is positive definite. \(G\) and \(R\) do not depend on \(c\).

The random-effect contribution is

G\Omega G^\mathsf{T} = \begin{pmatrix}0.85&0.50\\0.50&0.40\end{pmatrix}, \qquad V= \begin{pmatrix}0.95&0.50\\0.50&0.60\end{pmatrix}.

The covariance derivative is \(\Omega_c=\begin{pmatrix}0&1\\1&0\end{pmatrix}\). Only the middle product-rule term remains, so

\frac{\partial V}{\partial c} =G\Omega_cG^\mathsf{T} = \begin{pmatrix}1&1\\1&0\end{pmatrix} \begin{pmatrix}0&1\\1&0\end{pmatrix} \begin{pmatrix}1&1\\1&0\end{pmatrix} = \begin{pmatrix}2&1\\1&0\end{pmatrix}.

This derivative is symmetric. Its zero in the lower-right entry does not imply that the corresponding variance is zero; it means that this entry of \(V\) does not change to first order as \(c\) changes in this particular design and model.

Parameters Shared Across Components

Often a single parameter affects more than one part of \(V\). For instance, a model parameter can change the sensitivity matrix and also determine a random-effect variance or residual-error scale. The total derivative must then include every nonzero term. The following example makes all three routes visible.

Worked Example: Sensitivity, Random-Effect, and Residual Changes Together

Let \(\psi=1\) be dimensionless, \(G(\psi)=(\psi,\,1)^\mathsf{T}\), scalar \(\Omega(\psi)=\psi\), and \(R(\psi)=\psi^2\operatorname{diag}(1,2)\). At \(\psi=1\), \(\Omega=1>0\) and \(R=\operatorname{diag}(1,2)\), so both covariance components are valid. Compute the covariance and its derivative.

At \(\psi=1\), \(G=(1,1)^\mathsf{T}\), \(G_\psi=(1,0)^\mathsf{T}\), \(\Omega_\psi=1\), and \(R_\psi=2\operatorname{diag}(1,2)=\operatorname{diag}(2,4)\). The covariance is

V = \begin{pmatrix}1\\1\end{pmatrix} \begin{pmatrix}1&1\end{pmatrix} + \begin{pmatrix}1&0\\0&2\end{pmatrix} = \begin{pmatrix}2&1\\1&3\end{pmatrix}.

Substitution into the full derivative formula gives the sensitivity contributions, random-effect contribution, and residual contribution, respectively:

\begin{aligned} G_\psi\Omega G^\mathsf{T}+G\Omega G_\psi^\mathsf{T} &= \begin{pmatrix}2&1\\1&0\end{pmatrix},\\ G\Omega_\psi G^\mathsf{T} &= \begin{pmatrix}1&1\\1&1\end{pmatrix},\\ R_\psi&= \begin{pmatrix}2&0\\0&4\end{pmatrix}. \end{aligned}

Adding them gives \(\partial V/\partial\psi=\begin{pmatrix}5&2\\2&5\end{pmatrix}\). Each contribution is symmetric in this example, and their sum is symmetric. Leaving out the change in either \(\Omega\) or \(R\) would produce a different, incomplete derivative.

Parameterizations, Checks, and Failure Modes

The derivative is with respect to the parameter as it is defined in the model, not necessarily a variance or standard deviation on its natural scale. If a variance is parameterized as \(\omega=\exp(\lambda)\), then the chain rule gives \(\partial\Omega/\partial\lambda=(\partial\Omega/\partial\omega)(\partial\omega/\partial\lambda)\), with \(\partial\omega/\partial\lambda=\exp(\lambda)=\omega\). A derivative with respect to \(\lambda\) therefore differs from a derivative with respect to \(\omega\), even though both describe the same model through different coordinates.

A practical diagnostic is to compare the analytic derivative with a central finite difference at the same fixed design:

\frac{\partial V}{\partial\psi} \approx \frac{V(\psi+h)-V(\psi-h)}{2h}.

The approximation is useful for checking calculations, not a substitute for considering parameter constraints. Choose \(h\) small enough to limit truncation error but large enough to avoid severe rounding error, and ensure both perturbed parameter values define valid covariance components. Agreement should be checked across a range of reasonable step sizes. Automatic differentiation can also calculate derivatives through the model code, but it does not correct an incorrectly specified covariance expression or an omitted dependency.

Common warning signs include a nonsymmetric derivative matrix, an unexpected zero contribution, or a numerical derivative that changes substantially with small changes in step size. Inspect the dependency of each component on the parameter, verify dimensions in each matrix product, and check the parameterization and constraints. A derivative matrix need not itself be positive semidefinite: only the covariance at valid parameter values must satisfy the covariance constraints. For example, increasing a covariance parameter can increase some entries of \(V\) while decreasing others.

Approach	Strength	Limitation
Analytic product and chain rules	Shows each source of change and can be evaluated accurately.	Can be algebraically demanding when \(G\) depends on nonlinear model sensitivities.
Finite differences	Useful as an independent numerical check.	Step size, rounding, or invalid perturbed covariance parameters can distort the result.
Automatic differentiation	Can propagate derivatives through implemented model calculations.	Depends on correct code and parameter dependencies; it does not establish that the model expression is appropriate.

A common mistake is to differentiate only \(\Omega\), on the grounds that the parameter is called a random-effect variance. That reasoning is wrong whenever the same parameter also changes \(G\) or \(R\). The correction is to differentiate the complete expression for \(V\), identify every parameter-dependent factor, and apply the product and chain rules to each one. The resulting covariance derivatives are ingredients for later population information calculations, but they are not themselves a measure of information.

Check Your Understanding

Try each question, then compare your work with the solutions below.

If only \(R\) depends on \(\psi\), what is \(\partial V/\partial\psi\)?
Suppose \(G\) and \(R\) are fixed and \(\Omega_\psi=I\). What expression gives \(\partial V/\partial\psi\)?
Why does a symmetric covariance \(V\) have a symmetric derivative with respect to a scalar parameter, when it is differentiable?
Does a negative eigenvalue of \(\partial V/\partial\psi\) necessarily indicate an invalid covariance model? Explain.

Solutions. (1) The derivative is \(R_\psi\), because the sensitivity and random-effect covariance terms are constant. (2) It is \(GIG^\mathsf{T}\), from the middle product-rule term. (3) Each entry satisfies \(V_{ij}(\psi)=V_{ji}(\psi)\) for all valid \(\psi\); differentiating that equality gives \((\partial V/\partial\psi)_{ij}=(\partial V/\partial\psi)_{ji}\). (4) No. A derivative describes local change, not a covariance matrix; it need not be positive semidefinite. Validity is assessed for \(V\), \(\Omega\), and \(R\) at the parameter value, subject to the model’s covariance assumptions.