From a Mixed-Effects Model to a Computable Matrix
A population Fisher matrix is only as reliable as the derivatives used to build it. In a nonlinear mixed-effects model, both the population mean and the within-subject covariance can depend on the parameters. A numerical implementation must therefore differentiate both, assemble their contributions correctly, and check that the resulting matrix behaves as expected.
The calculation in this tutorial uses a Gaussian approximation to each subject's marginal outcome distribution. It is a local information calculation at specified nominal parameter values and a fixed design; it is not a guarantee that the approximation describes the exact marginal distribution of a nonlinear mixed-effects model.
This expression applies when \(m_i\) and \(V_i\) are differentiable in the parameters, and \(V_i\) is symmetric positive definite at the evaluation point. The design is treated as fixed. Subjects are independent of one another; measurements within a subject need not be independent, since their dependence is represented by \(V_i\). If the design has covariate strata, compute the appropriate subject-level contributions at each stratum and sum over the planned subjects.
Why Both Derivative Terms Appear
For one subject, set \(W=V^{-1}\) and \(z=Y-m\). Ignoring a constant, the Gaussian log likelihood is $$ \ell=-\frac{1}{2}\log|V|-\frac{1}{2}z^{T}Wz. $$ Differentiating with respect to \(\theta_a\) gives a mean part and a covariance part: $$ \frac{\partial\ell}{\partial\theta_a} =m_{,a}^{T}Wz+\frac{1}{2}\left(z^{T}WV_{,a}Wz-\operatorname{tr}(WV_{,a})\right). $$ Under the model, \(E(z)=0\) and \(E(zz^T)=V\). The expected product of two mean parts is \(m_{,a}^{T}Wm_{,b}\). The expected product of two centered covariance parts is \(\frac{1}{2}\operatorname{tr}(WV_{,a}WV_{,b})\), using the fourth-moment identity for a multivariate normal vector. The cross-products have expectation zero because they involve centered Gaussian third moments. Adding these terms gives the stated Fisher formula. In particular, the trace term is required whenever covariance depends on a parameter.
Worked Example: Analytic Derivatives for a Nonlinear Mixed Model
Worked Example: Analytic Derivatives for a Nonlinear Mixed Model
Consider a hypothetical model with two measurements per subject, at normalized times with \(g=(1,2)^T\). The conditional response is \(Y_{ij}=\exp(\beta+\eta_i)g_j+\epsilon_{ij}\), where \(\eta_i\sim N(0,q)\), the residual errors are independent with variance \(r\), and \(\eta_i\) is independent of the residual errors. Linearizing the response in \(\eta_i\) at zero gives the Gaussian approximation \(m=e^\beta g\) and \(V=qmm^T+rI\). Take \(\theta=(\beta,q,r)\) and evaluate at \((0,0.25,0.5)\). These values and the design are hypothetical.
At this point, \(m=(1,2)^T\) and $$ V=0.25 \begin{pmatrix}1&2\\2&4\end{pmatrix} +0.5 \begin{pmatrix}1&0\\0&1\end{pmatrix} = \begin{pmatrix}0.75&0.5\\0.5&1.5\end{pmatrix}. $$ Its determinant is \(0.75(1.5)-0.5^2=0.875\), so it is positive and $$ W=V^{-1}=\frac{1}{0.875} \begin{pmatrix}1.5&-0.5\\-0.5&0.75\end{pmatrix} =\frac{1}{7} \begin{pmatrix}12&-4\\-4&6\end{pmatrix}. $$
Because \(m_{,\beta}=m\), the derivatives are \(V_{,\beta}=2qmm^T\), \(V_{,q}=mm^T\), and \(V_{,r}=I\). At the evaluation point, the derivative matrices are $$ V_{,\beta}=\begin{pmatrix}0.5&1\\1&2\end{pmatrix},\qquad V_{,q}=\begin{pmatrix}1&2\\2&4\end{pmatrix},\qquad V_{,r}=\begin{pmatrix}1&0\\0&1\end{pmatrix}. $$ The mean term contributes only to \(I_{\beta\beta}\). Since \(Wm=(4/7)m\), it is \(m^TWm=20/7\).
For the covariance terms, \(V_{,\beta}=0.5V_{,q}\). Also, with \(A=V_{,q}=mm^T\), \(\operatorname{tr}(WAWA)=(m^TWm)^2=400/49\). Thus the \(q,q\) covariance contribution is \(200/49\), and the \(\beta,q\) and \(\beta,\beta\) covariance contributions are \(100/49\) and \(50/49\). For terms involving \(r\), \(m^TW^2m=320/49\), giving contributions \(160/49\) for \(q,r\) and \(80/49\) for \(\beta,r\). Finally, \(\frac12\operatorname{tr}(W^2)=106/49\) for \(r,r\). Combining the mean and covariance terms yields $$ I_i= \begin{pmatrix} 190/49&100/49&80/49\\ 100/49&200/49&160/49\\ 80/49&160/49&106/49 \end{pmatrix} \approx \begin{pmatrix} 3.878&2.041&1.633\\ 2.041&4.082&3.265\\ 1.633&3.265&2.163 \end{pmatrix}. $$ The entries are information on the specified parameter scale; their units follow from the units of the parameters. The matrix is symmetric, as required.
Computing and Checking Derivatives Numerically
Analytic derivatives are often efficient and accurate, but can be difficult to derive or maintain for a complicated model. A central finite difference provides a direct numerical check. For any scalar or matrix-valued function \(F\), approximate its derivative with respect to \(\theta_a\) by $$ F_{,a}(\theta)\approx \frac{F(\theta+h_a e_a)-F(\theta-h_a e_a)}{2h_a}, $$ where \(e_a\) is the unit vector for parameter \(a\), and \(h_a\) is a small perturbation in that parameter's units. Use the same parameter ordering for the mean and covariance derivatives as for the final matrix.
A practical implementation evaluates \(m_i\) and \(V_i\) at the nominal parameter vector, calculates or approximates every derivative, computes \(V_i^{-1}\), and fills each entry using the formula. It then adds subject contributions. Avoid explicitly forming matrix inverses in production code when a stable linear solve can be used: for example, solve systems involving \(V_i\) to evaluate the required products and traces.
Worked Example: Checking a Mean Derivative by Central Difference
Use the same hypothetical model, with \(g=(1,2)^T\) and nominal \(\beta=0\). Since \(m(\beta)=e^\beta g\), the analytic derivative is \(m_{,\beta}=g=(1,2)^T\). Set \(h=0.0001\). The central difference is $$ \frac{m(0.0001)-m(-0.0001)}{0.0002} =\frac{e^{0.0001}-e^{-0.0001}}{0.0002}g =\frac{\sinh(0.0001)}{0.0001}g \approx 1.0000000017g. $$
The numerical result is approximately \((1.000000002,2.000000003)^T\), close to the analytic derivative. This check tests the derivative implementation, not the validity of the Gaussian approximation itself.
Worked Example: Checking the Covariance Derivative
In the same model, \(V(\beta)=q e^{2\beta}gg^T+rI\), so the analytic derivative at \(\beta=0\) is \(V_{,\beta}=2qgg^T=0.5gg^T\). With \(q=0.25\) and \(h=0.0001\), the central difference is $$ \frac{V(h)-V(-h)}{2h} =q\frac{\sinh(2h)}{h}gg^T \approx 0.5000000033 \begin{pmatrix}1&2\\2&4\end{pmatrix}. $$
The analytic coefficient is \(0.5\); the central-difference coefficient differs by about \(3.3\times10^{-9}\). In particular, the off-diagonal derivative is approximately \(1.0000000067\) and the lower-right derivative is approximately \(2.000000013\), versus analytic values \(1\) and \(2\). This checks that the implementation differentiates the parameter-dependent covariance, not just the mean.
Replication and Implementation Checks
Worked Example: Summing Information Across Replicated Subjects
Suppose a hypothetical cohort contains three independent subjects with the same two-time design and the same nominal parameters as in the first example. Independence means their Fisher contributions add. Therefore, the cohort matrix is \(3I_i\), or $$ I_{\mathrm{cohort}}\approx \begin{pmatrix} 11.63&6.122&4.898\\ 6.122&12.24&9.796\\ 4.898&9.796&6.490 \end{pmatrix}. $$
For example, the \(\beta,\beta\) entry is \(3(190/49)=570/49\approx11.63\), and the \(q,r\) entry is \(3(160/49)=480/49\approx9.796\). A program that returns the single-subject matrix for this cohort has omitted the subject sum; one that returns three times this matrix has counted each subject three times.
Several checks are useful before interpreting a computed matrix:
- Symmetry: \(I_{ab}\) and \(I_{ba}\) should agree up to numerical rounding. A substantial discrepancy often indicates inconsistent derivatives or matrix multiplication order.
- Positive semidefiniteness: expected Fisher information is positive semidefinite. A materially negative eigenvalue can indicate a coding error, an invalid covariance matrix, or inaccurate numerical derivatives. Small negative eigenvalues near numerical precision can occur in computation.
- Finite-difference stability: repeat derivative checks with more than one sensible step size. A step that is too large has truncation error; one that is too small can suffer cancellation or solver noise. Parameters with different scales may require different \(h_a\).
- Replication: identical independent subjects must multiply the subject-level matrix by their number. This is a straightforward check of the aggregation code.
- Covariance validity: check that each \(V_i\) is symmetric positive definite at the nominal parameters before inversion or solving.
Choices, Limitations, and a Common Mistake
Analytic derivatives can be fast and precise when available, and they make the parameter dependence explicit. Central finite differences are simple to implement and valuable for verification, but their accuracy depends on step size and numerical stability. Automatic differentiation is another way to obtain derivatives from code while reducing finite-difference error, although it still requires checking the model implementation and parameterization.
The formula is exact for a Gaussian model whose mean and covariance are \(m_i(\theta)\) and \(V_i(\theta)\). For a nonlinear mixed-effects model, those quantities may instead come from an approximation, such as the first-order linearization used here. The computed matrix is then information under that Gaussian approximation. It can be sensitive to the nominal parameter values, the approximation, and parameter scaling. A singular or nearly singular matrix can indicate weak local identification, but does not by itself identify the cause.
Common mistake: retaining only \(m_{,a}^{T}V^{-1}m_{,b}\) because it resembles a familiar mean-model information formula. That omits information carried by changes in the covariance. In the example, \(q\) and \(r\) have zero mean derivatives, so the mean term alone would assign them no information at all. The correction is to compute \(V_{,a}\) for every parameter and include the trace term whenever those derivatives are nonzero.
Check Your Understanding
Answer each question, then compare with the solutions below.
- Which two derivative contributions appear in the Gaussian Fisher formula, and what model feature gives rise to the trace contribution?
- If \(N\) independent subjects have identical designs and identical parameter values, how does their total Fisher matrix relate to the single-subject matrix?
- Why is a central finite difference with an extremely small step not necessarily more accurate?
- In the worked nonlinear model, why would a mean-only information calculation fail to quantify information about \(q\)?
Solutions. (1) The contributions are \(m_{,a}^{T}V^{-1}m_{,b}\) and \(\frac12\operatorname{tr}(V^{-1}V_{,a}V^{-1}V_{,b})\). The trace term arises because the covariance changes with parameters. (2) Independence makes information additive, so \(I_{\mathrm{total}}=NI_i\). (3) Very small perturbations can make the difference \(F(\theta+h e_a)-F(\theta-h e_a)\) vulnerable to floating-point cancellation or solver noise; step-size stability should be checked. (4) The mean does not depend on \(q\), so its derivative with respect to \(q\) is zero. Information about \(q\) in this approximation comes through the covariance derivative \(V_{,q}\) and its trace contribution.