Why Numerical Validation Matters
Local influence summarizes how a fitted model responds to small changes under a specified perturbation scheme. Its curvature is calculated from derivatives of the perturbed log-likelihood at the unperturbed fit. A coding error, an inaccurate optimizer, or a poorly scaled perturbation can therefore produce a plausible-looking curvature result that does not describe the model’s actual behavior. Numerical checks ask whether the derivatives are being computed correctly and whether the local approximation predicts nearby refits.
The checks here apply to selection-model analyses, including settings with outcome-dependent dropout. They validate a calculation; they do not establish that the selection model is scientifically correct or that its missing-data assumptions are justified.
These definitions presume an interior fitted maximum, a smooth likelihood and perturbation scheme, and a locally identifiable parameterization. For the usual interpretation of curvature as nonnegative, \(H\) should be negative definite at the fitted maximum. In computation, solve linear systems involving \(H\) rather than explicitly forming \(H^{-1}\); direct inversion can amplify numerical error.
The curvature matrix follows from the first-order response of the fitted parameter. Differentiating the likelihood score equation \(\ell_\theta(\hat{\theta}(\omega);\omega)=0\) at the baseline gives \(H\,d\hat{\theta}+\Delta\,d\omega=0\). Thus \(d\hat{\theta}/d\omega=-H^{-1}\Delta\). Substituting this response into the second-order change in the baseline log-likelihood yields the quadratic approximation to likelihood displacement:
This is a local approximation, not an identity for finite changes. The numerical checks below evaluate both its ingredients and its predictive performance.
Finite-Difference Checks
Finite differences approximate derivatives by evaluating the log-likelihood at nearby parameter values. For a scalar coordinate with step size \(h\), the central difference for a first derivative is \(\{\ell(x+h)-\ell(x-h)\}/(2h)\). For a second derivative it is \(\{\ell(x+h)-2\ell(x)+\ell(x-h)\}/h^2\). For distinct coordinates \(i\) and \(j\), a common mixed-derivative approximation is:
The same four-corner formula can check entries of \(H\) or of \(\Delta\), using parameter coordinates or perturbation coordinates as appropriate. Derivative checks should use the same likelihood definition and parameterization as the fitted analysis. In particular, every evaluation must apply the same treatment of missing observations, constraints, and likelihood contributions.
Worked Example: Checking a Mixed Derivative
Suppose a hypothetical analysis of 40 patients with incomplete follow-up has a two-parameter log-likelihood contribution near its fitted value described by \(\ell(\theta_1,\theta_2)=-\tfrac12\theta_1^2-\theta_2^2+0.3\theta_1\theta_2+0.05\theta_1^2\theta_2\). The coordinates are centered at the fitted parameter, \((0,0)\). The mixed Hessian entry should be \(H_{12}=0.3\). Check it using equal steps \(h_1=h_2=0.1\).
The four-corner numerator is \(\ell(0.1,0.1)-\ell(0.1,-0.1)-\ell(-0.1,0.1)+\ell(-0.1,-0.1)\). The quadratic terms contribute zero to this signed sum. The \(0.3\theta_1\theta_2\) term contributes \(4(0.3)(0.1)(0.1)=0.012\). The \(0.05\theta_1^2\theta_2\) term cancels because its value is unchanged when \(\theta_1\) changes sign and the signed sum subtracts those values.
Therefore the finite-difference estimate is \(0.012/\{4(0.1)(0.1)\}=0.3\), matching the analytic derivative. If a numerical implementation instead returned, for example, 0.48, the discrepancy would warrant checking the evaluation points, coordinate order, and denominator before interpreting any curvature based on this Hessian.
A Hessian of a smooth scalar log-likelihood is symmetric: \(\ell_{ij}=\ell_{ji}\). Compare independently calculated off-diagonal entries, preferably using a consistent stencil and function-evaluation routine. With identical four-corner evaluations, the two entries may be algebraically identical; a separate calculation can still expose swapped indices or inconsistent step sizes. Large asymmetry is a warning about implementation, convergence, or numerical precision—not a substantive feature of a correctly evaluated smooth likelihood.
Reduce the Perturbation Size and Compare Refits
A finite-difference result depends on step size. If the step is too large, truncation error makes the finite difference differ from the derivative at the fitted point. If it is too small, floating-point cancellation and optimizer or likelihood-evaluation noise can dominate. There is no universally best step: assess stability over a sequence of decreasing steps, scaled to the magnitudes and units of the coordinates. For example, compare \(h\), \(h/2\), and \(h/4\). A sequence that settles to similar values is more reassuring than a single calculation.
Derivative stability alone does not establish that the local quadratic approximation describes a useful range of perturbations. For selected directions, refit the model at actual small perturbations. Compare the resulting likelihood displacement and parameter change with their local predictions. Use the same fitting criteria and check that each refit has converged to the intended solution. A change of optimizer, a different maximum, or a nearly singular Hessian can invalidate a superficial comparison.
Worked Example: Does the Curvature Approximation Improve as the Step Shrinks?
For a hypothetical fitted selection model with one scalar perturbation \(\omega\), suppose the profile log-likelihood is \(\ell(\theta;\omega)=-\tfrac12\theta^2-\tfrac14\theta^4+\omega\theta\), with baseline \(\omega_0=0\). The baseline estimate is \(\hat{\theta}_0=0\), \(H=-1\), and \(\Delta=1\), so \(K=-2(1)(-1)^{-1}(1)=2\). The local prediction is \(LD(\omega)\approx\omega^2\).
At a perturbed fit, the score equation is \(-\theta-\theta^3+\omega=0\), or \(\theta+\theta^3=\omega\). For \(\omega=0.1\), the solution is approximately \(\hat{\theta}(0.1)=0.09902\). The baseline log-likelihood at this refit is \(-\tfrac12\theta^2-\tfrac14\theta^4\), so \(LD=2\{0-\ell(\hat{\theta}(0.1);0)\}=\theta^2+\tfrac12\theta^4\approx0.009853\). The curvature prediction is \(0.1^2=0.01000\), close to the refit displacement.
For \(\omega=0.05\), solving the same equation gives \(\hat{\theta}(0.05)\approx0.04988\). Then \(LD\approx(0.04988)^2+\tfrac12(0.04988)^4\approx0.002491\), while the local prediction is \(0.05^2=0.002500\). The agreement improves as the perturbation shrinks, as expected for a valid local approximation. The curvature estimate calculated from a refit, \(2LD/\omega^2\), is about 1.971 at \(\omega=0.1\) and 1.993 at \(\omega=0.05\), approaching the derivative-based value \(K=2\).
Check the Curvature Matrix and Its Directional Predictions
When the perturbation has multiple coordinates, check that the calculated matrix \(K\) is symmetric and that its leading directions are stable to reasonable changes in finite-difference step. Under the stated assumptions, \(K\) is positive semidefinite: for any direction \(d\), \(d^TKd=-2(\Delta d)^TH^{-1}(\Delta d)\geq0\), because \(H^{-1}\) is negative definite. Tiny negative eigenvalues can arise from rounding. Materially negative eigenvalues suggest a failed assumption or a numerical problem, such as an indefinite Hessian, an inaccurate fit, or a coding error.
For a selected unit direction \(d\), compare \(LD(\omega_0+td)\) from an exact refit with \(\tfrac12t^2d^TKd\) for small positive and negative \(t\), when both are allowed by the perturbation scheme. A noticeable mismatch that decreases as \(|t|\) decreases is consistent with higher-order terms. A mismatch that persists or grows calls for investigation. Refits in a few directions are targeted diagnostics, not proof that the approximation is accurate everywhere.
Worked Example: Interpreting a Directional Refitting Check
Suppose a hypothetical longitudinal study gives a calculated curvature matrix \(K=\begin{pmatrix}8&2\\2&4\end{pmatrix}\) for two scaled perturbation coordinates. For the unit direction \(d=(1,0)^T\), the predicted curvature is \(d^TKd=8\), so the local displacement at step \(t\) is \(\tfrac12t^2(8)\). At \(t=0.05\), this equals \(0.0100\).
Suppose the exact refit at that perturbation gives \(LD=0.0106\). The relative difference from the local prediction is \((0.0106-0.0100)/0.0100=0.0600\), or 6.00%. At \(t=0.025\), the predicted displacement is \(\tfrac12(0.025)^2(8)=0.002500\); if the refit gives \(0.00254\), the relative difference is 1.60%. This shrinking discrepancy is consistent with a local approximation. The comparison is meaningful only if both refits converged and use the same likelihood and perturbation definition as the derivative calculation.
Diagnostics, Alternatives, and Limitations
A practical validation sequence is:
Check score residuals, optimizer convergence, and whether the fitted Hessian is negative definite and adequately conditioned.
Compare finite-difference derivatives across scaled step sizes and inspect Hessian symmetry.
Inspect symmetry, eigenvalues, and stability of leading directions under modest changes in numerical settings.
Compare exact likelihood displacement and parameter changes with local predictions at decreasing step sizes.
These checks complement rather than replace one another. Finite differences can reveal an incorrect derivative implementation but may themselves be unstable. Hessian symmetry is necessary for a smooth likelihood but does not guarantee accurate derivatives. Exact refits test local predictions at chosen points but can miss problems in untested directions. Analytic derivatives or automatic differentiation are useful alternatives to finite differences, especially when likelihood evaluations are noisy or the parameter dimension is large; they still require independent checks of the fitted model, derivative implementation, and refit behavior.
Common failure signals include a Hessian that changes sharply with step size, large asymmetry, a non-negative-definite Hessian at a claimed maximum, unstable leading curvature directions, refits that fail to converge, or discrepancies that do not diminish for smaller perturbations. Diagnose the source before reporting a ranking of influential subjects or perturbations. Near-boundary parameters, weak identification, flat likelihoods, and multiple local maxima can make all these calculations delicate. Finally, numerical agreement establishes only that the calculation is internally consistent near the fitted model; it does not validate the selection mechanism or resolve sensitivity to untestable missing-data assumptions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- If a central-difference Hessian has \(H_{12}=0.42\) and \(H_{21}=-0.10\), what should you investigate before calculating curvature?
- Why can making the finite-difference step extremely small make a derivative estimate worse rather than better?
- A unit direction has curvature \(C_d=6\). What is the local likelihood-displacement prediction at step \(t=0.02\)?
- Suppose an exact-refit discrepancy shrinks when the perturbation is halved. What does this support, and what does it not establish?
Solutions. (1) A smooth likelihood has a symmetric Hessian, so check coordinate ordering, stencils and step sizes, likelihood evaluation, and numerical precision; do not proceed as though the discrepancy were substantive. (2) Very small steps can cause cancellation when nearby likelihood values are subtracted, and evaluation or optimization noise can dominate the derivative signal. (3) The quadratic prediction is \(\tfrac12t^2C_d=\tfrac12(0.02)^2(6)=0.00120\). (4) It supports the interpretation that higher-order terms explain some of the discrepancy and that the local approximation is improving nearby. It does not prove accuracy in other directions or validate the missing-data model.