Gauss-Hermite Quadrature for Random-Effects Integration

Why integrate over random effects?

A dynamic prediction for an individual uses the longitudinal history available at a landmark time $t$ and the fact that the individual has remained event-free through $t$. The history updates what is plausible about that individual’s random effects. A prediction that averages over those updated possibilities is therefore generally different from one that simply substitutes a single estimated random-effects value.

The likelihood-based integration introduced in the earlier material on individual prediction leads to integrals over the random-effects distribution. Gauss–Hermite quadrature is a deterministic way to approximate these integrals. Its accuracy depends on using the right nodes and weights for the chosen normal distribution, and on checking that the quadrature adequately represents the integrand.

Definition: Let $b_i$ be an individual’s $q$-dimensional random-effects vector, with $b_i\sim N(0,D)$, where $D$ is positive definite. At fixed model parameters, let $L_i(t;b)$ denote the longitudinal-data likelihood from the individual’s observed marker history through $t$, conditional on random effects $b$. Let $S_i(v\mid b)$ denote conditional event-free survival through time $v$. For $u>t$, the landmark survival probability is $$ P(T_i>u\mid T_i>t,\mathcal Y_i(t)) = \frac{\int S_i(u\mid b)L_i(t;b)\phi_q(b;0,D)\,db} {\int S_i(t\mid b)L_i(t;b)\phi_q(b;0,D)\,db}, $$ where $\mathcal Y_i(t)$ is the observed longitudinal history and $\phi_q(\,\cdot\,;0,D)$ is the density of $N(0,D)$. The corresponding risk of an event in $(t,u]$ is one minus this probability.

This expression assumes that the joint model specifies the conditional marker and event contributions given $b$, that the random-effects distribution is correctly specified, and that the likelihood includes the observed marker history and survival to the landmark. It is conditional on the model parameters: integrating over $b$ does not, by itself, integrate over parameter uncertainty. The numerator is no greater than the denominator in the exact model because $S_i(u\mid b)\leq S_i(t\mid b)$ for $u>t$.

From a normal integral to quadrature

The standard-normal Gauss–Hermite rule is most simply written using nodes and weights for an expectation over $Z\sim N(0,1)$. If $D=LL^\mathsf{T}$ is a Cholesky factorization, then $b=LZ$ has distribution $N(0,D)$, where the components of $Z$ are independent standard normal variables. Thus a random-effects integral can be written as $$ \int f(b)\phi_q(b;0,D)\,db = E\{f(LZ)\}. $$

For one dimension, an $m$-node standard-normal quadrature rule has nodes $z_j$ and weights $w_j$, with weights summing to one. It approximates the expectation by $\sum_{j=1}^m w_j f(z_j)$. Many software references instead give the Gauss–Hermite rule for integrals weighted by $e^{-x^2}$. In that convention the standard-normal nodes are $z_j=\sqrt{2}x_j$, and the standard-normal weights are the $e^{-x^2}$ weights divided by $\sqrt{\pi}$. The scaling is essential.

Key result: For the three-node standard-normal rule, the nodes are $0,\sqrt{3},-\sqrt{3}$, with respective weights $2/3,1/6,1/6$. For $q$ independent standard-normal coordinates, the tensor-product rule uses every combination of the one-dimensional nodes, with a product weight. With random-effects covariance $D$, map each standard-normal vector $z$ to $b=Lz$ before evaluating the integrand.

The numerator and denominator can each be approximated by applying this rule to their integrands, then dividing the two sums. The same nodes and weights should be used in both sums. Since the ratio is computed from two approximations, it is important to assess the accuracy of both sums—not just whether the final prediction seems plausible.

Worked calculations

Worked Example: A Three-Node Landmark Prediction

Suppose a hypothetical individual is event-free at landmark $t$, and the prediction horizon is $u>t$. Consider a one-random-effect model and use three-node standard-normal quadrature. At the nodes $b=-\sqrt{3},0,\sqrt{3}$, suppose the hypothetical model evaluations of $g_t(b)=L_i(t;b)S_i(t\mid b)$ are $0.9588,0.9048,0.7884$, respectively. The evaluations of $g_u(b)=L_i(t;b)S_i(u\mid b)$ are $0.8814,0.7408,0.4901$. These are illustrative integrand values, not observations from a real study.

The denominator approximation uses weights $1/6,4/6,1/6$: $$ D_Q=\frac{0.9588+4(0.9048)+0.7884}{6} =\frac{5.3664}{6} =0.8944. $$ The numerator approximation is $$ N_Q=\frac{0.8814+4(0.7408)+0.4901}{6} =\frac{4.3347}{6} \approx 0.72245. $$

Therefore, the estimated conditional survival probability is $$ \widehat P(T_i>u\mid T_i>t,\mathcal Y_i(t)) =\frac{0.7224}{0.8944} \approx 0.8077, $$ and the estimated event risk in $(t,u]$ is $1-0.8077=0.1923$, or about $19.23\%$. The nodewise future integrand values are no greater than their landmark counterparts, and the computed numerator is below the denominator, as the survival ordering requires.

Worked Example: Comparing Quadrature Orders Against an Exact Value

Let $B\sim N(0,1)$ and consider $E\{\exp(-B^2/2)\}$. This is a simple diagnostic integral with an exact answer: $$ E\{\exp(-B^2/2)\} =\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-b^2}\,db =\frac{1}{\sqrt{2}} \approx 0.7071. $$

Using three standard-normal nodes, the center contributes $2/3$, and each of the two outer nodes has weight $1/6$. Since the function has value $e^{-3/2}$ at either outer node, $$ Q_3=\frac{2}{3}+\frac{1}{3}e^{-3/2} \approx 0.7410. $$ The absolute error is approximately $|0.7410-0.7071|=0.0339$.

For five standard-normal nodes, the center node is zero, the inner nodes are $\pm\sqrt{5-\sqrt{10}}$, and the outer nodes are $\pm\sqrt{5+\sqrt{10}}$. The center weight is $8/15$; each inner-node weight is $(7+2\sqrt{10})/60$; and each outer-node weight is $(7-2\sqrt{10})/60$. Substitution gives $$ Q_5=\frac{8}{15} +2\frac{7+2\sqrt{10}}{60}e^{-(5-\sqrt{10})/2} +2\frac{7-2\sqrt{10}}{60}e^{-(5+\sqrt{10})/2} \approx 0.7109. $$ The five-node absolute error is about $0.0038$, smaller than the three-node error in this example. This comparison illustrates a useful check, not a guarantee that error decreases monotonically for every integrand.

Worked Example: Mapping Correlated Random Effects

Suppose a hypothetical model has two random effects with covariance matrix $$ D=\begin{pmatrix}1&0.30\\0.30&0.25\end{pmatrix}, \qquad L=\begin{pmatrix}1&0\\0.30&0.40\end{pmatrix}, \qquad D=LL^\mathsf{T}. $$ For a tensor-product three-node rule, the standard-normal coordinates are each selected from $0,\pm\sqrt{3}$. For example, $z=(\sqrt{3},0)^\mathsf{T}$ maps to $$ b=Lz=(\sqrt{3},0.30\sqrt{3})^\mathsf{T} \approx (1.732,0.5196)^\mathsf{T}. $$ The node $z=(0,\sqrt{3})^\mathsf{T}$ maps to $b=(0,0.40\sqrt{3})^\mathsf{T}\approx(0,0.6928)^\mathsf{T}$. Each coordinate has three possible values, so there are $3^2=9$ tensor-product nodes; each node’s weight is the product of its coordinate weights.

As a check on the transformation and weights, consider $f(b)=\exp(a^\mathsf{T}b)$ with $a=(0.20,-0.40)^\mathsf{T}$. The exact expectation is $\exp(a^\mathsf{T}Da/2)$. Here $Da=(0.08,-0.04)^\mathsf{T}$, so $a^\mathsf{T}Da=0.032$, and the exact answer is $e^{0.016}\approx1.0161$.

Because $a^\mathsf{T}Lz=0.08z_1-0.16z_2$ and the tensor rule factorizes, the quadrature estimate is $$ \left\{\frac{2}{3}+\frac{1}{3}\cosh(0.08\sqrt{3})\right\} \left\{\frac{2}{3}+\frac{1}{3}\cosh(0.16\sqrt{3})\right\} \approx 1.0161. $$ This agrees closely with the exact value. The node mapping uses the full standard-normal coordinates $0,\pm\sqrt{3}$; omitting the standard-normal scaling would incorrectly shrink the random-effects nodes.

Implementation, accuracy, and failure modes

A reliable implementation keeps the statistical integral and the numerical approximation distinct. For each landmark and horizon, evaluate the history likelihood and survival contribution at each transformed node, form the two weighted sums, and then calculate the ratio. When likelihood contributions are very small, summing them directly can underflow. A practical implementation can work with log integrand values, combine each with the log quadrature weight, and use a log-sum-exp calculation before taking the ratio.

Set the target integral.
Write the landmark denominator and future numerator using the same observed history and random-effects distribution. Confirm that the numerator uses survival through $u$, not survival only through $t$.

Choose and transform the rule.
Use standard-normal nodes and weights, or convert the $e^{-x^2}$ convention by multiplying nodes by $\sqrt{2}$ and normalizing the weights. Map multivariate nodes using a Cholesky factor of $D$.

Evaluate both sums.
At each node, compute the longitudinal-history likelihood and the relevant conditional survival contribution. Apply the product weights in multiple dimensions.

Check numerical stability.
Repeat with a higher quadrature order, inspect numerator and denominator separately, and check that the resulting conditional survival probability lies between zero and one.

An important difficulty is that the integrand may be concentrated in a region that the fixed nodes do not represent well. This can happen when the marker history strongly updates the random-effects distribution, when the posterior mass is far from zero, or when the likelihood and survival terms create sharp curvature. Increasing the number of nodes is a basic diagnostic. If estimates shift materially as the order increases, the lower-order calculation is not reliable enough to report without further checks.

In higher dimensions, tensor-product rules become expensive: with $m$ nodes per coordinate and $q$ random effects, they require $m^q$ evaluations. A high-dimensional integral may therefore need a different strategy. Adaptive quadrature can place nodes more effectively around the important region, while Monte Carlo integration uses random draws and has simulation error. Laplace approximation is another alternative when a suitable local approximation is accurate. These approaches have different computational and approximation trade-offs; none removes the need to check the prediction.

Approach	Practical strength	Important limitation
Fixed Gauss–Hermite quadrature	Deterministic, with transparent nodes and weights; effective for smooth, low-dimensional normal integrals.	Tensor-product cost grows quickly with dimension, and fixed nodes can miss concentrated integrands.
Adaptive quadrature	Can focus evaluations around a shifted or concentrated region of integration.	Requires additional setup and can still be inaccurate if the integrand is poorly represented.
Monte Carlo integration	Can be more feasible as dimension grows.	Has simulation variability and requires enough draws to stabilize the risk estimate.
Laplace approximation	Can be computationally efficient when the integrand is well approximated locally.	Accuracy may suffer when the integrand is skewed, multimodal, or poorly captured by a local approximation.

A useful diagnostic is to compare predictions across several quadrature orders while holding the model and data fixed. Also check that the denominator is positive, that the numerator does not exceed it beyond rounding error, and that the predicted event risk remains in $[0,1]$. A violation can indicate incorrect weights, a node-scaling error, inconsistent integrands, or numerical instability. Agreement across orders is reassuring but not a proof of accuracy: different orders can share a systematic problem, so an exact test integral or an independent numerical calculation can provide a stronger check.

Common mistake and scope

A frequent error is to take nodes from the $e^{-x^2}$ convention, use them directly as standard-normal random effects, and then apply standard-normal weights. This mixes two distinct rules. In the three-node case, the unscaled nodes are $0,\pm\sqrt{3/2}$; the standard-normal nodes are $0,\pm\sqrt{3}$. The weights must be used in the matching convention as well. In multiple dimensions, each standard-normal node vector must also be multiplied by $L$ to represent covariance $D$.

Gauss–Hermite quadrature addresses numerical integration over random effects at fixed model parameters. It does not resolve misspecification of the longitudinal model, the event model, or the random-effects distribution. Nor does it automatically account for parameter uncertainty. The numerical approximation should be reported as part of the computational method, with enough information about quadrature order and accuracy checks for the prediction to be reproducible.

Key takeaway: For a normal random-effects integral, transform standard-normal quadrature nodes to the model’s random-effects scale, use the corresponding weights consistently in both landmark integrals, and verify accuracy by increasing the quadrature order or comparing against an independent check.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

For three-node standard-normal quadrature, what are the nodes and weights? What is the center weight?
If $D=LL^\mathsf{T}$, how should a standard-normal node vector $z$ be converted to a random-effects vector?
Why is it not enough to increase the quadrature order once and assume that the answer is accurate?
A quadrature calculation gives a numerator of $0.62$ and a denominator of $0.80$. What are the conditional survival probability and event risk?

Solutions. (1) The nodes are $0,\pm\sqrt{3}$; the weights are $2/3$ at zero and $1/6$ at each outer node, so the center weight is $2/3$. (2) Use $b=Lz$, which gives $b\sim N(0,D)$ when $z\sim N(0,I)$. (3) Convergence should be assessed by comparing several orders or an independent check; a single increase does not establish that the result has stabilized. (4) The conditional survival probability is $0.62/0.80=0.775$, and the event risk is $1-0.775=0.225$, or $22.5\%$ over the specified prediction interval.