Worked EM Iteration for a Shared-Random-Intercept Model

Why Work Through One EM Iteration?

In a joint model, each subject's random intercept influences both the repeated biomarker and the event hazard. The random intercept is unobserved, so an EM algorithm averages complete-data quantities over its conditional distribution given that subject's observed data. A single iteration makes this logic concrete: the E-step computes posterior moments, and the M-step uses them to update parameters.

The calculation below uses a deliberately small hypothetical data set and three-node Gaussian quadrature. It is a hand-checkable illustration, not a claim that three nodes are adequate for a real analysis. To keep the updates transparent, the association parameter is fixed at its current value during this iteration; in an application where it is estimated, its M-step must also maximize the expected complete-data log-likelihood.

Model, Data, and Assumptions

The random intercepts are independent across subjects. Censoring is assumed noninformative conditional on the modeled information, including \(b_i\), so the censoring mechanism contributes no additional parameter-dependent factor to the event-time likelihood. The event-time contribution, omitting factors that do not depend on the model parameters, is \( [\lambda\exp(\alpha b_i)]^{\delta_i}\exp\{-\lambda T_i\exp(\alpha b_i)\} \). The Gaussian model and the stated conditional-independence assumption are also part of the specification, not conclusions established by an EM calculation.

Use the following hypothetical observations and starting values. Subject 1 has \(Y_1=0\), an event at \(T_1=2\), and \(\delta_1=1\). Subject 2 has \(Y_2=0\), is censored at \(T_2=1.45\), and has \(\delta_2=0\). The starting values are \(\beta=0\), \(\sigma^2=1\), \(D=1\), \(\lambda=0.2\), and fixed \(\alpha=0.5\). Thus the longitudinal likelihood at a random-intercept value \(b\), with common constants removed, is \(\exp(-b^2/2)\) for both subjects.

The E-Step: Posterior Weights and Moments

For each subject, the posterior distribution of \(b_i\) is proportional to its random-effects density multiplied by its longitudinal and event-time likelihood contributions. With \(D=1\), a three-node quadrature rule for the standard normal uses nodes \(-\sqrt{3},0,\sqrt{3}\) and weights \(1/6,2/3,1/6\). For any function \(g\), the quadrature approximation to its posterior expectation is the sum of \(g\) at each node times that node's normalized posterior weight.

Here \(w_k\) are the prior quadrature weights, and \(L_i(b_k)\) is the subject's likelihood contribution at node \(b_k\), with factors common to all nodes omitted. The normalization is subject-specific: the two subjects generally have different posterior weights even though they share the same random-effects distribution.

For subject 2 there is no event-density factor, because the observation is right-censored. The survival contribution is \(\exp\{-0.2(1.45)\exp(0.5b)\}=\exp\{-0.29\exp(0.5b)\}\). Multiplying this by the Gaussian longitudinal likelihood and the prior quadrature weights gives unnormalized products approximately \(0.03292,0.49884,0.01866\). Their sum is approximately \(0.55042\); the normalized weights are approximately \(0.0598,0.9063,0.0339\). Hence the moments required below are \(\widehat E(b_2\mid\text{data})\approx -0.0449\), \(\widehat E(b_2^2\mid\text{data})\approx 0.2811\), and \(\widehat E\{\exp(0.5b_2)\mid\text{data}\}\approx 1.0121\). The hat indicates quadrature approximation, not a new parameter estimate.

The M-Step: Updating Parameters

The EM M-step maximizes the expected complete-data log-likelihood, where the expectation is under the posterior distribution of the latent random intercepts computed at the starting values. The posterior moments are held fixed while maximizing. In this example, there is one biomarker measurement per subject and the association parameter \(\alpha=0.5\) is held fixed, so the relevant updates have simple forms.

The first three equations follow by differentiating the expected Gaussian log-likelihood and the normal random-effects log-likelihood with respect to \(\beta,\sigma^2,D\). For the last equation, the survival part of the expected complete-data log-likelihood that depends on \(\lambda\) is \((\sum_i\delta_i)\log\lambda-\lambda\sum_iT_i\widehat E\{\exp(\alpha b_i)\mid\text{data}\}\). Setting its derivative to zero yields the displayed rate update. This maximizer is positive when there is at least one event and the denominator is positive.

Why Posterior Moments Matter—and How to Check Them

Several checks help locate errors in a hand calculation. Each set of normalized weights must be nonnegative and sum to one. A posterior second moment must be at least the square of its posterior mean; for subject 1, \(0.2876>(0.0734)^2\). The updated variances \(D^{\mathrm{new}}\) and \((\sigma^2)^{\mathrm{new}}\) must be nonnegative, and the updated exponential rate must be positive. Recalculate node products before normalizing: a small error in a tail product changes both posterior moments and downstream parameter updates.

A particularly common mistake is to treat the central-node weight as part of the second moment, or to forget the square of the nonzero node. The correct calculation is \((\sqrt{3})^2(w_-+w_+)=3(w_-+w_+)\), not \(w_-+w_+\). Another mistake is to use the event likelihood for a censored observation: a censored subject contributes survival through \(T_i\), but not the hazard factor at \(T_i\).

Three-node quadrature is useful here because every operation can be inspected, but it may be inaccurate when the posterior is skewed, concentrated in a tail, or otherwise poorly represented by the selected nodes. More nodes or another numerical integration scheme can improve the E-step approximation. Monte Carlo EM approximates the same posterior expectations by simulation and introduces simulation variability; direct maximum likelihood instead optimizes an integrated likelihood rather than using these EM updates. These approaches do not remove the need to check the model, the numerical integration, or the behavior of the fitted likelihood. In particular, a single EM step does not establish convergence, identifiability, or adequacy of the Gaussian and exponential assumptions.