Why Differentiate the Quadrature Approximation?
Direct maximum likelihood for a joint model often requires integrating a subject-specific likelihood over latent random effects. When that integral is approximated by quadrature, an optimizer needs the gradient of the approximated log likelihood. Differentiating the approximation carefully matters: its gradient is not generally the average of the node scores with the original quadrature weights, and parameters can affect the random-effects values at the nodes as well as the conditional likelihood there.
This tutorial focuses on that calculation and on checking it with finite differences. The likelihood is assumed to factor across subjects conditional on their random effects, with independent subject-level random effects drawn from a specified distribution. The conditional longitudinal and event-time contributions must be differentiable in the parameters being optimized. For the Gaussian random-effects calculations below, the covariance matrix is positive definite so it has a differentiable Cholesky factor locally. These conditions allow the usual chain rule to be applied to the quadrature sum. The earlier material on marginal likelihood differentiation provides the corresponding integral-level perspective.
Differentiate the Weighted Sum
For Gaussian random effects, write \(b_i=\mu(\theta)+A(\theta)z\), where \(A(\theta)A(\theta)^\mathsf{T}=D(\theta)\), the random-effects covariance matrix, and \(z\) has a standard multivariate normal distribution. Using fixed quadrature nodes \(z_k\) and weights \(w_k\), with weights normalized to sum to one, gives
The random-effects density is accounted for by the standard normal expectation and its quadrature weights; it is not an additional factor to multiply into \(q_{ik}\). Differentiating the log of the finite sum gives
The normalized values \(\rho_{ik}\) are nonnegative and sum to one. They are the relative contributions of the quadrature nodes to the approximated subject likelihood. The derivative at a node includes both the direct effect of \(\theta\) on the conditional likelihood and any indirect effect through the transformed random effect:
The first term holds \(b_i\) fixed. The second is the chain-rule contribution from how the random-effects value at the node changes with \(\theta\). If a parameter changes neither the conditional likelihood at fixed \(b_i\) nor the transformation, its node derivative is zero. If quadrature nodes or weights themselves depend on \(\theta\), their derivatives must also be included; the standard-normal quadrature formulation above uses fixed nodes and weights.
A Short Derivation and a Stable Calculation
For one parameter \(\theta_j\), differentiate the finite sum first. Since the weights are fixed, \(\partial\widehat L_i/\partial\theta_j=\sum_k w_k\,\partial q_{ik}/\partial\theta_j\). Using \(\partial q_{ik}/\partial\theta_j=q_{ik}\,\partial\log q_{ik}/\partial\theta_j\), then dividing by \(\widehat L_i\), yields
Thus the weights in the score are the normalized likelihood contributions, not just the quadrature weights. For independent subjects, the full approximated log likelihood is \(\widehat\ell(\theta)=\sum_i\log\widehat L_i(\theta)\), and its gradient is the sum of the subject-specific gradients. In computation, evaluating \(\log\widehat L_i\) with a log-sum-exp calculation helps avoid underflow when the node likelihoods are small. The same normalization should be used when calculating the score.
Worked Example: A Fixed-Effect Score
Consider one hypothetical longitudinal observation \(y=1\) with conditional model \(Y\mid b\sim N(\beta+0.5b,1)\), and random intercept \(b\sim N(0,\sigma^2)\). Take \(\beta=0.2\), \(\sigma=0.5\), and use the two-node standard-normal quadrature rule \(z_+=1,z_-=-1\), \(w_+=w_-=0.5\). Terms constant in the parameters can be omitted from the likelihood kernel. The transformed random effects are \(b_+=0.5\) and \(b_-=-0.5\).
At each node let \(e=y-\beta-0.5b\), so the likelihood kernel is \(q=\exp(-e^2/2)\). The residuals are \(e_+=1-0.2-0.25=0.55\) and \(e_-=1-0.2+0.25=1.05\). Thus \(q_+=\exp(-0.55^2/2)\approx0.859632\) and \(q_-=\exp(-1.05^2/2)\approx0.576228\).
At fixed \(b\), the derivative of \(\log q\) with respect to \(\beta\) is \(e\). The normalized node contributions are \(\rho_+=0.859632/(0.859632+0.576228)\approx0.5984\) and \(\rho_-\approx0.4016\). Therefore,
The quadrature-approximated score is about \(0.7507\) per unit of \(\beta\). It is positive at this parameter value: increasing the conditional mean slightly increases the approximated likelihood locally.
Worked Example: A Random-Effects Scale Score
Use the same hypothetical model, observations, parameter values, and nodes. Now differentiate with respect to \(\sigma\). Because \(b=\sigma z\), the node value changes with the scale: \(\partial b/\partial\sigma=z\). The conditional log-likelihood derivative with respect to \(b\) is \(0.5e\), not \(e\), since the conditional mean contains \(0.5b\). The full node score is consequently \(0.5ez\).
At the positive node, this score is \(0.5(0.55)(1)=0.275\). At the negative node it is \(0.5(1.05)(-1)=-0.525\). Using the same likelihood kernels as above gives
The approximated score is slightly negative, so the likelihood decreases locally as \(\sigma\) increases from \(0.5\). The sign and magnitude depend on the chain rule and on the relative node contributions. In particular, using \(e z\) for the scale score would omit the factor \(0.5\) in the conditional mean and give the wrong gradient.
Checking With Finite Differences
A finite-difference check compares the analytic gradient with numerical changes in the same quadrature-approximated objective. For a scalar parameter \(\theta_j\), the central difference is
where \(e_j\) is the unit vector for coordinate \(j\). For a smooth function, central differences have truncation error of order \(h^2\), while floating-point cancellation becomes more important as \(h\) gets very small. A useful check therefore tries several step sizes on a scale appropriate to the parameter, rather than trusting one arbitrary step. For a positive scale parameter, one can also check in its log-scale parameterization, remembering that the derivative with respect to \(\log\sigma\) equals \(\sigma\,\partial\widehat\ell/\partial\sigma\).
Worked Example: A Central-Difference Check
For the one-observation model above, the approximated likelihood is \(\widehat L(\beta,\sigma)=\{q_+(\beta,\sigma)+q_-(\beta,\sigma)\}/2\), and the constant factor \(1/2\) cancels in the log-likelihood difference. At \(\beta=0.2,\sigma=0.5\), the log of the sum of the two node kernels is approximately \(0.3618\). Keep \(\beta=0.2\) fixed and check the scale score with \(h=0.001\).
At \(\sigma=0.501\), the residuals are \(1-0.2-0.5(0.501)=0.5495\) and \(1-0.2+0.5(0.501)=1.0505\). At \(\sigma=0.499\), they are \(0.5505\) and \(1.0495\). Evaluating the same two-node sum at these values gives log sums of approximately \(0.361719\) and \(0.361811\), respectively. Hence
The displayed log sums are rounded too coarsely for this small difference to reproduce the analytic score reliably. Using unrounded evaluations gives a central difference of about \(-0.0461\), matching the analytic scale score to the reported precision. This illustrates an important practical point: retain full numerical precision during the check, and round only when reporting. A mismatch that remains with full precision across sensible step sizes calls for checking the derivative, parameter transformation, and quadrature implementation.
Practical Choices, Diagnostics, and Limitations
Analytic differentiation, automatic differentiation, and finite differences serve different roles. Analytic differentiation makes the chain rule explicit and can be efficient, but it is vulnerable to algebraic omissions. Automatic differentiation applies the chain rule to the implemented finite sum and can reduce hand-calculation errors, but it differentiates the code as written; a coding error or an incorrect likelihood remains an error. Finite differences are useful as an independent check, but they require repeated evaluations and depend on a suitable step size.
| Approach | Strength | Important caution |
|---|---|---|
| Analytic gradient | Efficient and interpretable | Check every direct and transformation derivative |
| Automatic differentiation | Systematically applies the chain rule to the computation | Does not correct a misspecified or incorrectly coded objective |
| Finite differences | Independent numerical check requiring only objective evaluations | Results can be distorted by step size and rounding |
Several diagnostics help isolate problems. First, compare analytic and central-difference gradients for multiple parameters at an interior parameter value, using more than one step size. Second, verify that each node score includes the derivative of \(b_i=\mu+A z_k\) whenever the mean or covariance parameters change. Third, check that the normalized node weights sum to one and that the score is their weighted average. When using a log scale for a standard deviation, make sure the optimizer's gradient is with respect to the log scale rather than the original scale.
A gradient check does not establish that the quadrature approximation is accurate. It verifies differentiation of that approximation. An accurate derivative of a coarse approximation can still differ from the derivative of the exact integrated likelihood. Conversely, a small discrepancy between finite differences and the analytic score can result from numerical cancellation rather than a faulty formula. The finite-difference step should be large enough to avoid cancellation but small enough to reflect local changes; checking a sequence of decreasing steps reveals whether the result stabilizes.
Check Your Understanding
Answer the questions, then compare your work with the brief solutions below.
- For a fixed-node quadrature approximation, what weights combine the node-specific log-likelihood scores into the gradient of the log of the sum?
- If \(b=\sigma z\), what additional derivative is needed for a parameter \(\sigma\) that affects the conditional likelihood only through \(b\)?
- Why should a finite-difference check evaluate the same quadrature approximation used to compute the analytic gradient?
- For the worked scale example, why is the positive-node scale score \(0.275\), rather than \(0.55\)?
- What does agreement between analytic and finite-difference gradients establish, and what does it not establish?
Solutions. (1) Use the normalized contributions \(\rho_{ik}=w_kq_{ik}/\sum_\ell w_\ell q_{i\ell}\). (2) Include \((\partial\log q/\partial b)(\partial b/\partial\sigma)\), with \(\partial b/\partial\sigma=z\). (3) Otherwise the check compares derivatives of different numerical objectives, so disagreement may reflect different approximations rather than a gradient error. (4) The conditional mean contains \(0.5b\), so \(\partial\log q/\partial b=0.5e\); at the positive node \(e=0.55,z=1\), giving \(0.5(0.55)(1)=0.275\). (5) It supports correct differentiation of the implemented approximation at the checked point; it does not show that the quadrature approximation is close to the exact integrated likelihood.