Why Gauss–Hermite Quadrature Appears in Joint Models
In a joint model, a subject’s longitudinal measurements and event outcome are conditionally modeled given that subject’s random effect. To obtain the marginal likelihood contribution, the conditional contribution must be averaged over the random-effect distribution. Even with a one-dimensional random effect, this integral will generally have no closed form because the conditional likelihood can involve nonlinear mean functions, survival terms, or both.
Gauss–Hermite quadrature approximates this particular kind of integral by evaluating its integrand at a finite set of carefully chosen nodes and combining those values with weights. The essential practical issue is the normal-density scaling: standard Gauss–Hermite nodes are designed for the weight function \(e^{-x^2}\), not directly for a normal density. A change of variables connects the two.
In many shared-random-effects models, \(b_i\) is a random intercept with \(\mu=0\) and variance \(\sigma^2\). The formula below also allows a nonzero mean, which is useful for understanding the transformation. It applies to a single Gaussian random effect; it does not by itself handle a vector of correlated effects.
From a Normal Integral to Gauss–Hermite Form
The \(K\)-point Gauss–Hermite rule for the physicists’ weight function approximates an integral of the form \(\int_{-\infty}^{\infty} e^{-x^2}f(x)\,dx\) by \(\sum_{j=1}^K w_j f(x_j)\). Here \(x_j\) are the Gauss–Hermite nodes and \(w_j\) are their corresponding weights. The weights are positive and sum to \(\sqrt{\pi}\). They are weights for integration against \(e^{-x^2}\), not probability weights that sum to one.
For the normal random effect, set \(b=\mu+\sqrt{2}\sigma x\). Then \(db=\sqrt{2}\sigma\,dx\), and the normal density becomes \(\phi(b;\mu,\sigma^2)=\{1/(\sqrt{2\pi}\sigma)\}\exp[-(b-\mu)^2/(2\sigma^2)]\). Substitution gives
The factor \(1/\sqrt{\pi}\) normalizes the weights for an expectation. The transformed evaluation points are \(\mu+\sqrt{2}\sigma x_j\); changing the variance changes their spread, while changing the mean shifts all of them. For the common zero-mean case, this reduces to evaluations at \(\sqrt{2}\sigma x_j\).
A useful justification for the method is that the \(K\)-point Gauss–Hermite rule integrates \(e^{-x^2}f(x)\) exactly whenever \(f\) is a polynomial of degree at most \(2K-1\). This follows from the Gaussian quadrature construction: the nodes are roots of the degree-\(K\) orthogonal polynomial for the weight \(e^{-x^2}\), and the associated weights make the rule exact through degree \(2K-1\). This exactness does not mean that a general joint-model likelihood is integrated exactly; its transformed integrand need not be a polynomial of that degree.
Worked Calculations
Worked Example: A Three-Point Approximation to a Conditional Contribution
Suppose \(b_i\sim N(0,0.2^2)\), and a hypothetical subject’s conditional contribution is \(g_i(b)=\exp(-1.2+0.4b-0.5b^2)\). Use the three-point rule to approximate \(L_i\).
The transformed nodes are \(0\) and \(\pm\sqrt{3}(0.2)\), so they are approximately \(0\), \(0.3464\), and \(-0.3464\). The normalized weights at these points are \(2/3\), \(1/6\), and \(1/6\), respectively.
At the center, \(g_i(0)=e^{-1.2}\approx0.3012\). At the positive node, the exponent is \(-1.2+0.4(0.3464)-0.5(0.3464)^2\approx-1.1214\), so \(g_i(0.3464)\approx0.3258\). At the negative node, the exponent is \(-1.2-0.4(0.3464)-0.5(0.3464)^2\approx-1.3986\), so \(g_i(-0.3464)\approx0.2470\).
Therefore, \(L_i\approx(1/6)(0.2470)+(2/3)(0.3012)+(1/6)(0.3258)\). The three terms are approximately \(0.04117\), \(0.2008\), and \(0.05430\), giving \(L_i\approx0.2963\). This is an approximation to the marginal contribution, not the conditional contribution evaluated at \(b=0\).
Worked Example: Averaging a Nonlinear Event Contribution
For a hypothetical subject, suppose \(b_i\sim N(0,0.4^2)\), and the conditional event contribution after accounting for the observed follow-up is \(g_i(b)=\exp\{-0.8\exp(0.5b)\}\). Use the two-point Gauss–Hermite rule to approximate the marginal contribution. This example isolates the numerical integration step; \(g_i\) stands for a specified conditional contribution.
For two points, the nodes are \(-1/\sqrt{2}\) and \(1/\sqrt{2}\), and both weights are \(\sqrt{\pi}/2\). After normal-density scaling, the evaluation points are \(\sqrt{2}(0.4)(\pm1/\sqrt{2})=\pm0.4\). The normalized weights are therefore \(1/2\) at each point.
At \(b=-0.4\), \(g_i(-0.4)=\exp\{-0.8e^{-0.2}\}=\exp(-0.6550)\approx0.5194\). At \(b=0.4\), \(g_i(0.4)=\exp\{-0.8e^{0.2}\}=\exp(-0.9771)\approx0.3762\). Thus \(L_i\approx(0.5194+0.3762)/2=0.4478\). The contribution is an average under the random-effect distribution; evaluating only at its mean would give \(g_i(0)=e^{-0.8}\approx0.4493\), a different approximation.
Worked Example: Scaling When the Random-Effect Mean Is Nonzero
Let \(B\sim N(0.3,0.5^2)\), and take \(g(b)=\exp(-b^2/2)\). Use the three-point rule and compare it with the analytic expectation, which is available here as a check.
The three evaluation points are \(0.3\) and \(0.3\pm\sqrt{3}(0.5)\), or approximately \(-0.5660\), \(0.3000\), and \(1.1660\). The function values are approximately \(g(-0.5660)=0.8520\), \(g(0.3000)=0.9560\), and \(g(1.1660)=0.5067\). Applying the normalized three-point weights gives \(L\approx(1/6)(0.8520)+(2/3)(0.9560)+(1/6)(0.5067)=0.8638\).
Completing the square in the normal integral gives the exact expectation \(E\{\exp(-B^2/2)\}=(1+0.5^2)^{-1/2}\exp\{-0.3^2/[2(1+0.5^2)]\}\approx0.8628\). The three-point result is close in this example, but a small difference is expected because the transformed function is not a polynomial of degree five or less.
Using the Approximation in a Joint-Model Likelihood
For a one-dimensional shared random effect, the subject-specific marginal likelihood contribution has the form \(L_i=\int g_i(b)\phi(b;\mu,\sigma^2)\,db\), where \(g_i(b)\) is the full conditional contribution for that subject. It includes the relevant longitudinal and event-time factors conditional on the same random effect. Apply quadrature to this integral separately for each subject, then form the sample log likelihood as \(\ell=\sum_i\log(L_i)\).
The order of these operations matters. The subject’s likelihood is an integral of the conditional contribution, so the quadrature sum approximates \(L_i\); the logarithm is taken after that sum. In general, \(\log\{\sum_j a_j\}\) is not equal to \(\sum_j\log(a_j)\). Taking logarithms of node-specific contributions and adding them would therefore not approximate the intended marginal log likelihood.
When \(g_i(b)\) is a product of several conditional factors, evaluate the product at each transformed node before applying the weights. In computation, products of very small factors may underflow. A stable implementation can compute the weighted sum on the log scale using the log-sum-exp identity, while retaining the same nodes, weights, and normal-density scaling.
Diagnostics, Alternatives, and Limitations
Quadrature order is a practical accuracy check. Recompute the likelihood with increasing numbers of nodes and compare the subject-level contributions, total log likelihood, and fitted quantities that matter for the analysis. Agreement across several adequate orders supports numerical stability; a large change suggests that the integral is not yet well approximated. The integral may be especially challenging if the conditional contribution changes sharply over regions with appreciable random-effect density, or if it is concentrated far from the normal mean.
Other useful checks include confirming that the normal distribution’s standard deviation—not its variance—is used in the node transformation, checking that normalized weights sum to one, and verifying that the result is finite and nonnegative when the integrand is a likelihood. For a positive likelihood contribution, the quadrature result should also be positive because the quadrature weights are positive. If the computed value is zero, infinite, or unstable, check both scaling and numerical underflow or overflow.
| Approach | Strength | Limitation |
|---|---|---|
| Gauss–Hermite quadrature | Deterministic and often efficient for a smooth, one-dimensional Gaussian integral. | Accuracy depends on node count and how well the fixed nodes cover the important part of the integrand. |
| Monte Carlo integration | Uses random draws and can be adapted to integrals that are difficult for a fixed grid. | Has simulation error, so estimates can vary between runs and may require many draws for precision. |
| Analytic integration | Can give an exact marginal contribution when the model’s form permits it. | Often unavailable for the nonlinear conditional contributions encountered in joint models. |
Fixed-node Gauss–Hermite quadrature is not guaranteed to be accurate merely because the integrand is smooth: the nodes and weights are concentrated according to the chosen normal reference distribution. For sharply concentrated or shifted contributions, checking higher orders is important. Adaptively centering nodes around the posterior mode is a distinct strategy and is treated separately.
Common Mistakes
- Using the variance where the standard deviation belongs. The transformed nodes use \(\sqrt{2}\sigma x_j\), not \(\sqrt{2}\sigma^2 x_j\).
- Forgetting the normalization. The raw Gauss–Hermite weights sum to \(\sqrt{\pi}\). For the normal expectation, use \(w_j/\sqrt{\pi}\), or use the already normalized coefficients in the three-point formula.
- Using unscaled nodes. The standard nodes \(x_j\) are not the random-effect values. They must be transformed to \(\mu+\sqrt{2}\sigma x_j\).
- Taking the logarithm at each node and summing. First approximate the integral to obtain \(L_i\); then take \(\log(L_i)\).
- Assuming a low-order rule is exact for any likelihood. Polynomial exactness applies only to a specific degree in the transformed integrand, not to an arbitrary nonlinear contribution.
Check Your Understanding
Try each question, then compare your reasoning with the short solutions below.
- For \(b\sim N(0,0.3^2)\), what is the three-point quadrature location corresponding to the positive nonzero standard Gauss–Hermite node?
- Why do the normalized three-point weights sum to one even though the raw Gauss–Hermite weights sum to \(\sqrt{\pi}\)?
- For a positive conditional likelihood contribution, what does a computed negative quadrature approximation indicate?
- In what order should the quadrature sum and logarithm be applied to obtain a subject’s marginal log-likelihood contribution?
Solutions. (1) The location is \(0+\sqrt{3}(0.3)=0.5196\), approximately. (2) The raw weights are divided by \(\sqrt{\pi}\); equivalently, the three normalized weights are \(1/6,2/3,1/6\), which sum to one. (3) It indicates a computational or implementation error, since the quadrature weights and conditional likelihood values are positive. (4) First compute the weighted sum approximating \(L_i\), then compute \(\log(L_i)\).