Constructing the EM Auxiliary Function

Why the Auxiliary Function Must Be Constructed Carefully

In a shared-random-effects joint model, each subject’s longitudinal measurements and event history are conditionally related through latent random effects. The observed-data likelihood integrates those effects out, whereas the complete-data likelihood treats them as if they were observed. The EM algorithm uses the latter to form an objective that can be optimized more easily, but its expectation must be taken using the random-effects posterior at the current parameter values.

This tutorial focuses on constructing that expected complete-data log likelihood, usually called the EM auxiliary function. It is not obtained by simply plugging a posterior mean into the complete-data likelihood: nonlinear terms, including survival contributions, generally require expectations over the full posterior distribution.

Complete-Data Likelihood and the E-Step Distribution

For subject \(i\), let \(y_{ij}\) be the longitudinal outcome at time \(t_{ij}\), for \(j=1,\ldots,n_i\). Let \(T_i\) be the observed event or censoring time and let \(\delta_i=1\) indicate an event and \(\delta_i=0\) indicate right censoring. Let \(b_i\) be the subject’s vector of shared random effects, and write \(O_i\) for all observed data on that subject, including the longitudinal and event-time data and the covariates.

Write \(\theta\) for the full parameter vector. A standard shared-random-effects specification assumes conditional independence of the longitudinal observations and event process given \(b_i\) and covariates. Its subject-specific complete-data log likelihood is

Here \(f_Y\) is the longitudinal conditional density, \(h_i\) is the event hazard, \(H_i(t\mid b_i;\theta)=\int_0^t h_i(u\mid b_i;\theta)\,du\) is the cumulative hazard, and \(f_b\) is the random-effects density. The survival expression is the log event-time contribution for either an observed event or right censoring; when \(\delta_i=0\), the hazard-at-event term is absent. This form assumes censoring is noninformative for the event-time model conditional on the covariates and random effects, or otherwise appropriately handled in the likelihood.

At an EM iteration, denote the current parameter value by \(\theta^{\mathrm{old}}\). The E-step distribution for subject \(i\) is

It is the posterior based on all of that subject’s observed data under the current parameters, not just the longitudinal measurements. The auxiliary function is the sum of expectations of the complete-data log likelihood under these distributions:

The candidate parameter \(\theta\) appears inside the complete-data log likelihood, while the distribution \(q_i\) stays fixed at the current iteration. In expanded form, the expectation applies separately to the longitudinal log-density terms, the survival terms, and the random-effects log-density. Terms constant in \(\theta\) may sometimes be dropped when optimizing, but retaining them makes the construction and comparisons transparent.

What the Expectation Does to Each Model Component

For a Gaussian longitudinal model with conditional mean \(m_{ij}(b_i;\theta)\) and residual standard deviation \(\sigma\), one contribution is

Thus the needed expectation is a posterior expected squared residual. In a linear random-intercept model, it can be calculated from the posterior mean and variance; in more complex models, it may require numerical integration. For survival, the expectation includes both \(E_{\mathrm{old}}[\delta_i\log h_i(T_i\mid b_i;\theta)\mid O_i]\) and \(E_{\mathrm{old}}[H_i(T_i\mid b_i;\theta)\mid O_i]\). The latter is not generally equal to the cumulative hazard evaluated at the posterior mean of \(b_i\).

Why This Construction Supports EM

The role of \(Q\) follows from a likelihood identity. Let \(p_\theta(b\mid O)\) be the posterior of the random effects under a candidate \(\theta\). For any fixed density \(q(b)\) supported where this posterior is positive, the Kullback–Leibler divergence is \(D_{\mathrm{KL}}(q\|p_\theta)=\int q(b)\log\{q(b)/p_\theta(b\mid O)\}\,db\), which is nonnegative. Since \(p_\theta(O,b)=p_\theta(O)p_\theta(b\mid O)\), rearranging gives

Set \(q\) equal to the posterior at \(\theta^{\mathrm{old}}\). The entropy term \(-\int q\log q\) is then constant as \(\theta\) varies. At \(\theta^{\mathrm{old}}\), the divergence is zero. Consequently,

Because the final term is nonnegative, increasing \(Q\) from its value at the current parameter cannot decrease the observed-data log likelihood. This justification depends on using the actual current posterior for the E-step and the same complete-data model in both the expectation and the observed-data likelihood. The identity explains why constructing \(Q\) correctly matters; it does not require that every expectation have a closed form.

Approximations, Diagnostics, and Common Pitfalls

When posterior expectations are not available analytically, they may be evaluated by numerical integration, such as quadrature over the random effects. The earlier material on Gauss–Hermite quadrature and posterior-focused integration provides the numerical integration tools; here the key requirement is to apply the approximation to the expectation of the complete-data log-likelihood terms under the current posterior. Check that posterior weights are normalized, that estimates are stable as integration accuracy is increased, and that each term uses the same posterior approximation.

Direct maximization of the integrated observed-data likelihood is a main alternative: it targets the same likelihood without constructing an EM auxiliary function, but it requires evaluating the integrated likelihood as parameters vary. EM instead separates posterior expectation from maximization. Its auxiliary function may still require numerical integration, and it need not have closed-form updates. The purpose of this tutorial is to construct \(Q\), not to derive parameter-specific M-step updates.

Common warning signs include a Q value that changes when only the candidate parameter is changed inside the E-step posterior, quadrature weights that do not sum to one, and a survival expectation computed only at the posterior mean despite a nonlinear hazard. A further diagnostic is to compare the component-wise sum with a direct numerical expectation of the full complete-data log likelihood at a test parameter.

• Wrong reasoning: “The random effect is unobserved, so I can replace it everywhere by its posterior mean and evaluate the complete-data log likelihood there.”
• Correction: The auxiliary function is an expectation of the log likelihood over the full current posterior. Posterior variance and nonlinear functions of the random effects can affect that expectation, as the worked calculations show.
• Wrong reasoning: “For each candidate parameter, I should recompute the posterior inside \(Q\).”
• Correction: In the E-step, form the posterior using \(\theta^{\mathrm{old}}\). Hold it fixed while evaluating \(Q(\theta\mid\theta^{\mathrm{old}})\); update the posterior at the next iteration.