The EM Algorithm With Latent Random Effects

Why Treat Random Effects as Latent Data?

In a joint model, each subject contributes longitudinal measurements and a time-to-event outcome, while also having an unobserved subject-specific random effect. The random effect links repeated measurements within a subject and can also link the longitudinal trajectory to the event hazard. Since it is unobserved, the likelihood for the observed data must integrate it out. That integration can be difficult, especially with correlated, multivariate random effects.

The expectation-maximization (EM) algorithm offers a different way to organize maximum-likelihood estimation. It treats the random effects as missing data: if their values were known, the likelihood would often have a simpler form. At each iteration, EM takes an expectation of that complete-data log likelihood under the current conditional distribution of the random effects, then optimizes the resulting expectation over the model parameters. This is a framework for estimation, not a promise that every update has a closed-form solution.

Definition: Let \(D_i=(y_i,T_i,\delta_i)\) denote subject \(i\)'s observed longitudinal data, event or censoring time, and event indicator, respectively. Let \(b_i\) be the subject-specific random effect and \(\theta\) the vector of model parameters. The observed-data likelihood integrates out \(b_i\); the complete-data likelihood treats \(b_i\) as if it were observed.

Likelihoods and Assumptions

Suppose subjects are independent, and the random effects have a specified density \(p(b_i;\theta)\). Write \(f(D_i\mid b_i;\theta)\) for the conditional density or likelihood contribution of the observed outcomes given the random effect. Then the subject's marginal likelihood contribution and the sample likelihood are

L_i(\theta)=\int f(D_i\mid b_i;\theta)p(b_i;\theta)\,db_i, \qquad L(\theta)=\prod_{i=1}^{n}L_i(\theta).

The integral is over the dimension of \(b_i\). In the complete-data formulation, the random effects are included alongside the observed data. The complete-data log likelihood is

\ell_c(\theta;D,b) = \sum_{i=1}^{n} \left\{ \log f(D_i\mid b_i;\theta)+\log p(b_i;\theta) \right\}.

This expression includes both the conditional outcome model and the random-effects distribution. Omitting the latter would discard the part of the likelihood that estimates parameters such as the random-effects covariance. In a common specification, conditional on \(b_i\) and covariates, the longitudinal and event-time contributions factor as \(f_Y(y_i\mid b_i;\theta)f_S(T_i,\delta_i\mid b_i;\theta)\). That factorization is a modeling assumption, not a universal property of joint models; use the conditional density specified by the model at hand.

For a time-to-event contribution, a standard form is \(h_i(T_i\mid b_i;\theta)^{\delta_i}S_i(T_i\mid b_i;\theta)\), where \(h_i\) is the hazard and \(S_i\) the survival function. This assumes the censoring mechanism is ignorable for the parameters being estimated, such as when censoring is conditionally independent and its parameters are distinct. The longitudinal and event models, random-effects distribution, and censoring assumptions must all be specified appropriately before applying EM.

The basic EM setup also requires a well-defined conditional distribution of the latent effects given the observed data and current parameters, and finite expectations of the complete-data log likelihood. For continuous random effects these expectations may themselves require numerical integration. The algorithm does not remove the integration problem; it changes the quantity that must be evaluated at each iteration.

The E-Step and M-Step

Let \(\theta^{(r)}\) be the parameter value at iteration \(r\). The conditional density of the random effects given the observed data under this value is

p(b\mid D;\theta^{(r)}) = \frac{f(D\mid b;\theta^{(r)})p(b;\theta^{(r)})} {L(\theta^{(r)})}.

The numerator combines the outcome likelihood conditional on \(b\) with the random-effects density. The denominator normalizes it to integrate to one. This conditional distribution is the basis of the E-step; its calculation for continuous random effects can be the numerically demanding part.

Definition: The EM auxiliary function at iteration \(r\) is the conditional expectation of the complete-data log likelihood, with the expectation taken under the current parameter value:

Q(\theta\mid\theta^{(r)}) = E_{\theta^{(r)}}\!\left[ \ell_c(\theta;D,b)\mid D \right].

In the E-step, form or evaluate this expectation as a function of the candidate parameter \(\theta\). In the M-step, choose \(\theta^{(r+1)}\) to maximize \(Q(\theta\mid\theta^{(r)})\), or at least increase it when exact maximization is impractical. The expectation uses the old parameter value to describe uncertainty about the latent effects, while the argument \(\theta\) in \(Q\) is the value being optimized. Keeping these two roles separate is essential.

A short derivation clarifies why the conditional distribution enters. By the definition of a conditional density, \(f(D,b;\theta^{(r)})=L(\theta^{(r)})p(b\mid D;\theta^{(r)})\). Thus, at the current fit, the latent effects are weighted according to how plausible each value is after observing both the longitudinal and event-time data. The E-step averages the complete-data log likelihood using those weights; it does not substitute a single estimated random effect for every subject.

Worked Example: Separating the Complete-Data Contributions

Consider one hypothetical subject with a scalar random intercept \(b\). Suppose \(b\sim N(0,0.36)\), the longitudinal observation is \(y=0.4\) with conditional model \(Y\mid b\sim N(0.2+b,0.25)\), and the event hazard is constant at \(h=\lambda\exp(\alpha b)\), with \(\lambda=0.2\) per year and \(\alpha=0.8\). The subject is censored at \(T=1.5\) years. For illustration, evaluate the complete-data log contribution at \(b=0.3\).

The random-effects log density is \(-\frac12\log(2\pi(0.36))-\frac{0.3^2}{2(0.36)} \approx -0.5331\). The longitudinal mean at this value of \(b\) is \(0.2+0.3=0.5\), so its log density contribution is \(-\frac12\log(2\pi(0.25))-\frac{(0.4-0.5)^2}{2(0.25)} \approx -0.2458\).

The hazard is \(0.2\exp(0.8(0.3))\approx0.2542498\) per year. For censoring at 1.5 years, the event-time likelihood contribution is the survival probability \(S(T\mid b)=\exp(-hT)\); its log is \(-hT=-(0.2\exp(0.24))(1.5)\approx-0.3814\). Therefore, at this fixed value of \(b\), the complete-data log contribution is approximately \(-0.5331-0.2458-0.3814=-1.1603\).

This is not the observed-data log likelihood: the random effect is not known to equal 0.3. The observed contribution integrates the joint density over all possible \(b\). In EM, values of \(b\) receive weights from their conditional distribution given this subject's observed data.

A Finite-Support Illustration of the E-Step

For continuous random effects, the conditional expectation may require quadrature or another integration method. A finite-support example makes the weighting logic visible without numerical integration. It is only a simplified illustration; standard joint models often use continuous random effects.

Worked Example: Conditional Weights for Two Latent Effects

Suppose a hypothetical model has a random effect \(b\) taking values \(+\) and \(-\), with \(P(b=+)=\pi=0.4\). Conditional on \(b=+\), a binary longitudinal measurement has probability \(0.8\) of being 1; conditional on \(b=-\), that probability is \(0.2\). The event hazard is constant at \(0.2\) per year for \(b=+\) and \(0.5\) per year for \(b=-\). Assume conditional independence of the measurement and event time given \(b\). The subject has measurement \(y=1\) and is censored at \(T=1\) year.

For censoring, the conditional event-time likelihood is the survival probability. The joint likelihoods conditional on the two possible random effects are therefore \(L_+=0.8\exp(-0.2(1))\approx0.6549846\) and \(L_-=0.2\exp(-0.5(1))\approx0.1213061\).

The unnormalized posterior weights are the prior probabilities multiplied by these likelihoods: \(w_+=0.4(0.6549846)\approx0.2619938\) and \(w_-=0.6(0.1213061)\approx0.0727837\). Their sum is approximately \(0.3347775\). Hence the conditional probability of \(b=+\) is \(0.2619938/0.3347775\approx0.7826\), and the probability of \(b=-\) is approximately \(0.2174\). These normalized weights—not a hard assignment to \(+\) or \(-\)—are used to average complete-data quantities in the E-step.

An Expected-Count M-Step

Worked Example: Updating a Random-Effects Probability

Continue the finite-support model, now with two hypothetical independent subjects. The first has the data just described: \(y_1=1\), censored at 1 year. The second has \(y_2=0\), also censored at 1 year. Keep the conditional measurement probabilities and hazards fixed, and estimate only \(\pi=P(b=+)\). Start from \(\pi^{(r)}=0.4\).

For subject 1, the E-step weight for \(b_1=+\) is \(\gamma_1\approx0.7826\), as calculated above. For subject 2, the likelihoods conditional on \(b=+\) and \(b=-\) are \(0.2\exp(-0.2)\approx0.1637462\) and \(0.8\exp(-0.5)\approx0.4852245\). The unnormalized weights are \(0.4(0.1637462)=0.0654985\) and \(0.6(0.4852245)=0.2911347\). Thus \(\gamma_2=P(b_2=+\mid D_2;\pi^{(r)})\approx0.1837\).

The part of the expected complete-data log likelihood involving \(\pi\) is \(\gamma_1\log\pi+\gamma_2\log\pi+(2-\gamma_1-\gamma_2)\log(1-\pi)\). Maximizing this expression gives the expected fraction of subjects in the \(+\) group: \[ \pi^{(r+1)}=\frac{\gamma_1+\gamma_2}{2} \approx\frac{0.7826+0.1837}{2}=0.4831. \] The update uses fractional expected membership, not a count based on assigning each subject to its most likely group. In a continuous random-effects model, the M-step instead involves expectations of functions of the effects and may require numerical optimization.

Diagnostics, Alternatives, and Limitations

EM is useful when integrating out random effects directly is difficult but expectations of complete-data quantities can be evaluated or approximated. It also provides a clear conceptual separation: the E-step accounts for uncertainty in the latent effects, and the M-step updates the model parameters using the resulting expected complete-data criterion. It does not, by itself, make either step computationally easy.

Monitor the observed-data log likelihood, parameter changes, and relevant numerical diagnostics across iterations. A fit may appear stable because parameter updates are tiny even when the likelihood calculation is inaccurate. For continuous effects, check that the numerical integration is sufficiently accurate and repeat the fit with tighter tolerances when feasible. Compare starting values, inspect estimates near parameter-space boundaries, and check that the fitted model respects covariance and hazard constraints. Poorly identified parameters, weak event information, or a nearly singular random-effects covariance can lead to slow progress or unstable estimates.

Exact EM with an exact E-step and a valid M-step is distinct from an approximate implementation. If the E-step uses inaccurate quadrature, simulation, or a poorly converged approximation, the computed update may not behave like the exact algorithm. Likewise, an M-step that fails to improve its objective can interrupt progress. When evaluating an implementation, distinguish numerical error from a model that simply fits poorly.

Approach	What it does	Practical trade-off
EM with latent random effects	Optimizes an expected complete-data log likelihood iteratively.	Can simplify parameter updates, but requires conditional expectations and may converge slowly near the maximum.
Direct marginal maximum likelihood	Optimizes the likelihood after integrating out the random effects.	Targets the marginal likelihood directly, but requires reliable integration and optimization together.
Laplace approximation	Approximates the random-effects integral near a mode.	Can be efficient, but accuracy depends on the quality of the approximation.
Quadrature-based likelihood	Approximates the integral using weighted evaluation points.	Can be accurate with suitable points, but computational cost can increase rapidly with dimension.

These approaches are alternatives for handling the same marginal likelihood, not different definitions of the target parameter. EM's main limitations are the cost of its E-step, potentially slow convergence, and the need to calculate or approximate expectations. It is also not a remedy for misspecification: incorrect conditional outcome models, an unsuitable random-effects distribution, or an inadequate censoring assumption can still yield misleading estimates.

Common mistake: using a best guess for each random effect. Replacing each \(b_i\) by its posterior mode ignores uncertainty and is not the EM E-step. The E-step averages the complete-data log likelihood over the full conditional distribution.
Common mistake: dropping the random-effects density. The complete-data likelihood includes \(\log p(b_i;\theta)\), which can contain parameters being estimated.
Common mistake: confusing the two parameter roles in \(Q\). The posterior distribution in the expectation is evaluated at \(\theta^{(r)}\); the candidate \(\theta\) is the argument optimized in the M-step.

Check Your Understanding

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

Why does the observed-data likelihood integrate over \(b_i\), while the complete-data likelihood includes its value?
In the finite-support example, why is the E-step weight for \(b=+\) not simply the prior probability \(0.4\)?
If three subjects have posterior probabilities of \(b_i=+\) equal to \(0.2\), \(0.6\), and \(0.9\), what is the M-step update for \(\pi\) when these are the only unknown parameter?
Name two assumptions or numerical features that should be checked when fitting a joint model by EM.

Solutions. (1) The random effects are unobserved, so the observed likelihood sums or integrates their joint density with the observed outcomes over possible values. The complete-data likelihood is defined as if those latent values were also observed. (2) The posterior weight reflects both the prior probability and how compatible each value is with the longitudinal and event-time data; here it is proportional to \(P(b=+)L_+\). (3) The update is \((0.2+0.6+0.9)/3=0.5667\), the average expected membership. (4) Examples include conditional independence assumptions in the specified model, the censoring assumption, accuracy of numerical integration in the E-step, and convergence from multiple starting values.