Numerical Integration Within the EM E-Step

Why Numerical Integration Enters the E-Step

In a joint model, each subject’s random effects are unobserved. The EM algorithm therefore forms an expectation over their distribution conditional on that subject’s observed longitudinal and event data, evaluated at the current parameter estimate. Even when the random effects have a Gaussian prior, the posterior density is generally not Gaussian: the survival contribution can make its shape asymmetric or concentrated. The required expectation usually has no closed form.

Numerical quadrature replaces the integral by a weighted sum. In the E-step, that sum must approximate a posterior expectation, not merely an expectation under the random-effects prior. The observed-data contribution changes the weights separately for each subject.

Definition: Let \(b_i\) be subject \(i\)’s random-effects vector, with prior density \(p(b_i;\theta)\), and let \(L_i(b_i;\theta)\) be that subject’s observed-data likelihood conditional on \(b_i\). At current parameter value \(\theta^{(r)}\), the E-step posterior density is \(p(b_i\mid \mathcal O_i;\theta^{(r)}) \propto L_i(b_i;\theta^{(r)})p(b_i;\theta^{(r)})\), where \(\mathcal O_i\) denotes the subject’s observed data. For any function \(g\) whose posterior expectation exists, \(\mathbb E[g(b_i)\mid\mathcal O_i;\theta^{(r)}]\) is the integral of \(g\) against this normalized posterior density.

Here \(L_i\) includes the longitudinal likelihood and the event-time likelihood, including the appropriate contribution for censoring. The likelihood and prior are evaluated at the current parameter values in the E-step. The resulting expectations supply terms in the EM auxiliary function, such as posterior means, second moments, and expectations involving the survival association structure.

From a Posterior Integral to a Weighted Sum

Suppose first that \(b_i\) is scalar and has a normal prior \(N(0,\sigma_b^2)\). A three-node Gauss–Hermite rule for an expectation under a standard normal distribution uses nodes \(x_k=(-\sqrt{3},0,\sqrt{3})\) and weights \(a_k=(1/6,2/3,1/6)\). For the stated variance, transform each node to \(b_k=\sigma_b x_k\). More generally, if \(b_i\sim N(0,D)\), take a matrix \(A\) satisfying \(AA^\mathsf T=D\), and transform a vector of standard-normal nodes \(x_k\) to \(b_k=Ax_k\).

For each node, multiply its prior quadrature weight by the subject’s conditional likelihood. Normalize these products to obtain approximate posterior weights. The posterior expectation of \(g\) is then approximated by the weighted average of \(g\) evaluated at the nodes.

\widehat{w}_{ik} = \frac{a_k L_i(b_k;\theta^{(r)})} {\sum_{\ell}a_\ell L_i(b_\ell;\theta^{(r)})}, \qquad \mathbb E[g(b_i)\mid\mathcal O_i;\theta^{(r)}] \approx \sum_k \widehat{w}_{ik}g(b_k).

The denominator approximates the marginal likelihood contribution for subject \(i\). In multiple dimensions, a tensor-product rule uses combinations of univariate nodes, with each combination’s prior weight equal to the product of its component weights. The number of nodes grows rapidly with dimension, which is one reason to choose an integration rule suited to the model’s random-effects dimension.

A Short Derivation of the Posterior-Weight Formula

The posterior expectation is a ratio of two integrals. For a scalar random effect, its numerator integrates \(g(b_i)L_i(b_i;\theta^{(r)})p(b_i;\theta^{(r)})\), and its denominator integrates \(L_i(b_i;\theta^{(r)})p(b_i;\theta^{(r)})\). Apply the same quadrature rule to both integrals. The common density transformation and integration constants cancel in the ratio; the remaining numerator and denominator are weighted sums.

\frac{\int g(b)L_i(b;\theta^{(r)})p(b;\theta^{(r)})\,db} {\int L_i(b;\theta^{(r)})p(b;\theta^{(r)})\,db} \approx \frac{\sum_k a_k L_i(b_k;\theta^{(r)})g(b_k)} {\sum_k a_k L_i(b_k;\theta^{(r)})} = \sum_k \widehat{w}_{ik}g(b_k).

This also explains why the likelihood belongs in the posterior weights. Using only the prior weights would calculate a prior expectation, which is not the E-step expectation conditional on the subject’s data.

Worked Example: Posterior Mean and Second Moment

Consider a hypothetical subject with scalar random effect \(b\sim N(0,1)\). Suppose the conditional observed-data likelihood is proportional to \(\exp\{-(1-b)^2/2\}\). Use the three-node rule above to approximate the posterior mean and second moment. The omitted proportionality constant is independent of \(b\), so it cancels when the weights are normalized.

The nodes are \(-\sqrt{3},0,\sqrt{3}\), approximately \(-1.7321,0,1.7321\). The likelihood values are approximately \(0.02394,0.6065,0.7649\). Multiplying by the prior quadrature weights gives unnormalized weights approximately \(0.0039906,0.40435,0.12749\), whose sum is \(0.53584\). Dividing by that sum gives posterior weights approximately \(0.0074475,0.7546230,0.2379296\).

The posterior mean is the weighted sum of the nodes:

(-1.7321)(0.0074475)+(0)(0.7546230)+(1.7321)(0.2379296) \approx 0.3992.

The second moment is the weighted sum of the squared nodes. Since the two nonzero nodes have square \(3\), this is \(3(0.0074475+0.2379296)\approx0.7361\). The posterior variance can be derived from these moments as approximately \(0.7361-(0.3992)^2=0.5767\). These are quadrature approximations to posterior quantities, not prior moments.

Expectations Needed Beyond the Posterior Mean

An E-step can require more than \(\mathbb E(b_i\mid\mathcal O_i)\). For a Gaussian random-effects covariance update, it typically needs the matrix second moment \(\mathbb E(b_i b_i^\mathsf T\mid\mathcal O_i)\). For a survival model with a current-value association, it can require a posterior expectation of a nonlinear function of the longitudinal trajectory. The same normalized posterior weights are used for each quantity: replace \(g\) by the function needed in the auxiliary function.

Worked Example: A Nonlinear Posterior Expectation

Continue with the hypothetical posterior weights and nodes in the preceding example. Suppose a survival-model term requires the posterior expectation of the random-effect multiplier \(\exp(0.5b)\). At the nodes, this function is approximately \(0.4206,1,\) and \(2.3774\). The weighted sum is

(0.0074475)(0.4206)+(0.7546230)(1)+(0.2379296)(2.3774) \approx 1.3234.

Thus the quadrature approximation is about \(1.323\), a dimensionless posterior expectation of the multiplier. It is not the exponential of the posterior mean: \(\exp\{0.5\,\mathbb E(b\mid\mathcal O_i)\}\) is a different quantity. In general, nonlinear functions must be evaluated at each node before averaging.

Worked Example: Scaling Nodes for a Different Prior Variance

Suppose a hypothetical random intercept has prior \(b\sim N(0,0.25)\), and the E-step requires its prior expectation of \(\exp(0.7b)\), with no data contribution in this calculation. The standard deviation is \(0.5\), so the three transformed nodes are \(-0.5\sqrt{3},0,0.5\sqrt{3}\), or approximately \(-0.8660,0,0.8660\). Their function values are approximately \(0.5454,1,1.8335\).

Using weights \(1/6,2/3,1/6\), the quadrature result is

\frac{1}{6}(0.5454)+\frac{2}{3}(1)+\frac{1}{6}(1.8335) \approx 1.0632.

The value is dimensionless. This example illustrates the variance transformation: using standard-normal nodes without multiplying by \(0.5\) would integrate under the wrong prior. The rule provides an approximation for this non-polynomial function; it is not exact merely because three nodes were used.

Numerical Checks, Stability, and Alternatives

Quadrature is a deterministic approximation, and its reliability should be assessed rather than assumed. A practical diagnostic is to recompute the E-step quantities using a higher-resolution rule and compare the results. Check posterior means, second moments, nonlinear expectations, and the resulting auxiliary-function terms—not just the marginal likelihood. A likelihood approximation may appear stable while a tail-sensitive posterior expectation is not.

For a multivariate random effect, tensor-product quadrature can become expensive as the dimension rises. Sparse-grid quadrature, covered in the earlier material on sparse-grid quadrature, can reduce the node count while retaining a structured deterministic rule. If the posterior is sharply concentrated away from where prior-centered nodes provide adequate coverage, increase resolution or consider whether prior-centered quadrature is appropriate for the problem.

Computations should also be protected against underflow. If likelihood values are extremely small, calculate log contributions \(\log(a_k)+\log L_i(b_k;\theta^{(r)})\), subtract their maximum, exponentiate, and normalize. This rescaling does not change normalized posterior weights. Verify that the normalized weights are nonnegative and sum to one up to rounding error.

Laplace approximation is a principal alternative: it uses a local quadratic approximation near a posterior mode. It can be efficient when that approximation captures the posterior shape, but can be less reliable for skewed or heavier-tailed posteriors. Quadrature instead samples a fixed integration rule and can directly approximate the required expectations, at a computational cost that increases with dimension and resolution. Simulation-based integration is another broad alternative; it replaces a deterministic weighted sum with a sample-based estimate and introduces simulation variability. The choice depends on dimension, posterior shape, computational resources, and the accuracy required for the E-step.

Common Mistakes and Limitations

Using prior weights as posterior weights. Prior weights alone ignore the subject’s observations. Multiply each prior quadrature weight by that subject’s conditional likelihood, then normalize.
Forgetting to transform nodes. For a nonstandard Gaussian prior, transform standard-normal nodes using the random-effects standard deviation or covariance factor.
Exponentiating an averaged random effect. For a nonlinear function \(g\), compute \(g(b_k)\) at each node and then average. In general, \(\mathbb E[g(b_i)\mid\mathcal O_i]\ne g\{\mathbb E(b_i\mid\mathcal O_i)\}\).
Assuming more nodes guarantee accuracy. More nodes often improve an approximation, but convergence should be checked for the posterior quantities actually used. A poorly chosen integration rule can remain inadequate at a modest resolution.

Quadrature also has limits. It approximates integrals; it does not remove model misspecification, and finite-resolution results are not exact E-step expectations. Its cost can be substantial for high-dimensional random effects, and a subject whose likelihood is concentrated in a region poorly represented by the chosen nodes can be especially challenging. These limitations motivate explicit integration checks and careful interpretation of the computed E-step.

Key takeaway: For each subject, form posterior quadrature weights by multiplying prior node weights by the conditional observed-data likelihood and normalizing. Use those weights for every required posterior moment, transform nodes to the correct covariance scale, and check stability as the integration resolution changes.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

Why must the conditional likelihood enter the quadrature weights in an E-step?
If \(b\sim N(0,4)\), how should standard-normal quadrature nodes be transformed?
Why is the posterior expectation of \(\exp(b)\) generally not the exponential of the posterior mean?
Name two quantities, besides the posterior mean, whose quadrature stability could be checked.

Solutions. (1) The E-step conditions on the observed data, so its expectation uses the posterior distribution, proportional to the prior density times the conditional likelihood. (2) The standard deviation is \(2\), so multiply each standard-normal node by \(2\). (3) The exponential is nonlinear: the function must be evaluated at each posterior node and then averaged, which generally differs from applying the function to the average. (4) Examples include the posterior second moment, the posterior variance, and the expectation of a nonlinear survival-model function.