Monte Carlo EM for Joint Models

Why Use Monte Carlo in the E-Step?

In a joint longitudinal and time-to-event model, a subject’s random effects are latent. The E-step averages the complete-data log likelihood over their posterior distribution, conditional on that subject’s observed longitudinal measurements and event or censoring information. Numerical integration can be difficult when the random-effects vector has several components or the posterior is sharply concentrated. Monte Carlo EM (MCEM) replaces the required posterior integrals with averages over simulated random effects.

This replacement makes the E-step practical in settings where direct numerical integration is burdensome, but it introduces simulation error. That error can affect the estimated auxiliary function, the M-step maximizer, and the apparent progress from one iteration to the next. The key is to distinguish error from the Monte Carlo approximation from the behavior of exact EM.

From the Posterior E-Step to a Monte Carlo Average

For subject $i$, write $D_i=(Y_i,T_i,\delta_i)$, where $Y_i$ denotes the observed longitudinal data, $T_i$ the observed follow-up time, and $\delta_i$ the event indicator. Let $b_i$ be the subject’s latent random-effects vector and $\theta$ the model parameters. The complete-data log likelihood contribution, including the random-effects distribution, is $\ell_{c,i}(\theta;D_i,b_i)$. Given the current parameter value $\theta^{(k)}$, the E-step uses the posterior distribution of $b_i$ given $D_i$ and $\theta^{(k)}$.

Definition: The EM auxiliary function at iteration $k$ is $$ Q(\theta\mid\theta^{(k)}) = \sum_{i=1}^{n} E_{\theta^{(k)}}\!\left[ \ell_{c,i}(\theta;D_i,b_i)\mid D_i \right]. $$ In MCEM, this expectation is approximated using draws from the posterior distribution $p(b_i\mid D_i;\theta^{(k)})$. With $M_i$ draws $b_i^{(1)},\ldots,b_i^{(M_i)}$, the approximation is $$ \widehat Q(\theta\mid\theta^{(k)}) = \sum_{i=1}^{n}\frac{1}{M_i} \sum_{m=1}^{M_i} \ell_{c,i}(\theta;D_i,b_i^{(m)}). $$

The draws are from the posterior at the current parameter value, not from the random-effects distribution before seeing the subject’s data. The longitudinal observations and survival information both contribute to that posterior in a joint model. In the M-step, maximize the estimated function with respect to $\theta$, holding the simulated draws fixed while carrying out that maximization. At the next EM iteration, the posterior distribution changes with the updated parameter value, so new draws are generated for the next E-step.

The basic MCEM description assumes that, conditional on the observed data and current parameters, the draws are independent samples from the required posterior distributions. The posterior expectations must exist, and the sampled quantities used to estimate them must have finite variance for the usual Monte Carlo standard-error calculation to apply. Subjects are treated as independent contributions under the model’s usual independent-subject sampling assumption. If the simulation method produces dependent draws, the independent-draw standard-error formula is not appropriate without adjustment; dependence reduces the effective amount of simulation information.

Simulation Error and Its Standard Error

For one subject, suppose the E-step requires $E[g(b_i)\mid D_i;\theta^{(k)}]$, where $g$ is a scalar function of the random effects. From $M$ independent posterior draws, estimate this expectation by the sample mean $\bar g$. If the sampled values have sample standard deviation $s_g$, the estimated Monte Carlo standard error (MCSE) is $s_g/\sqrt{M}$. This measures uncertainty due to using a finite simulation sample; it is not the standard error of a model parameter based on the observed study sample.

Worked Example: Estimating a Posterior Mean and Its MCSE

Suppose four independent hypothetical draws from one subject’s current random-effects posterior are $1.25,\ 2.25,\ 2.75,\ 3.75$. Estimate the posterior mean of the random effect and its MCSE.

The estimated posterior mean is $$ \bar b=\frac{1.25+2.25+2.75+3.75}{4} =\frac{10.00}{4}=2.50. $$ The deviations from the mean are $-1.25,\ -0.25,\ 0.25,\ 1.25$. Their squared sum is $$ (-1.25)^2+(-0.25)^2+(0.25)^2+(1.25)^2 =1.5625+0.0625+0.0625+1.5625=3.25. $$ Thus the sample standard deviation and estimated MCSE are $$ s_b=\sqrt{\frac{3.25}{4-1}}=\sqrt{1.0833}\approx1.0408, \qquad \operatorname{MCSE}(\bar b)=\frac{1.0408}{\sqrt{4}}\approx0.5204. $$ The estimate of the posterior mean is $2.50$ random-effects units, with an estimated simulation error of about $0.5204$ units. Four draws provide a quite imprecise estimate; this MCSE should not be confused with uncertainty about the underlying model parameters.

For an entire auxiliary function, uncertainty depends on the complete-data log-likelihood contributions being averaged. If contributions are independent across draws, the sample variance of those contributions gives an MCSE for the corresponding estimated average. For a parameter update that is a nonlinear function of several estimated expectations, an MCSE for each expectation alone does not automatically give an MCSE for the update. Replicate simulations or a suitable approximation can be used to assess that propagated uncertainty.

How Monte Carlo Error Enters an M-Step

The estimated auxiliary function can have a slightly different shape from the exact $Q$ function. As a result, its maximizer may differ from the exact EM update. This matters especially when $Q$ is relatively flat near its maximum: a small perturbation to the estimated function can move the estimated maximizer noticeably. A well-defined, sharply peaked maximum is generally less sensitive to the same perturbation.

Worked Example: A Monte Carlo M-Step for a Regression Parameter

Consider a simplified complete-data contribution in which an observed value $z=1$, conditional on one latent random effect $b$, has distribution $N(\alpha b,1)$. The parameter of interest is $\alpha$; other terms in the complete-data log likelihood do not depend on $\alpha$. Suppose four hypothetical posterior draws are $b=0,1,2,3$. Find the MCEM update for $\alpha$.

Ignoring terms constant in $\alpha$, the estimated auxiliary function is $$ \widehat Q(\alpha) =-\frac{1}{2(4)}\sum_{m=1}^{4}(1-\alpha b^{(m)})^2. $$ Differentiating and setting the derivative to zero gives $$ \frac{d\widehat Q}{d\alpha} =\frac{1}{4}\sum_{m=1}^{4}b^{(m)}(1-\alpha b^{(m)})=0, \qquad \widehat\alpha =\frac{\sum_m b^{(m)}}{\sum_m (b^{(m)})^2}. $$ For these draws, $\sum_m b^{(m)}=0+1+2+3=6$, and $\sum_m (b^{(m)})^2=0+1+4+9=14$. Therefore, $$ \widehat\alpha=\frac{6}{14}=\frac{3}{7}\approx0.4286. $$ This is the maximizer of the Monte Carlo approximation, not necessarily the maximizer of the exact auxiliary function. The discrepancy depends on how well these four draws represent the current posterior expectations of $b$ and $b^2$.

In this example, an unusual simulation batch could alter the ratio of the estimated moments and therefore the update. More generally, a finite simulation sample can produce a step that appears to reduce the observed-data likelihood, even though exact EM has a non-decrease property under its usual conditions. That property applies to an exact E-step followed by a valid M-step; it does not guarantee that a noisy estimate of $Q$ will give a non-decreasing likelihood at every iteration.

Choosing Simulation Effort and Diagnosing Trouble

For independent draws with a stable variance, the MCSE of a sample mean decreases at rate $M^{-1/2}$. Multiplying the number of draws by four therefore approximately halves its MCSE, although the realized sample standard deviation can also change. This slow rate makes simulation allocation important. Small batches may be adequate when far from convergence, but near convergence the likelihood gains and parameter changes can be smaller than the simulation noise.

Worked Example: Comparing Two Batch Sizes

Suppose an MCEM M-step for a location parameter $\gamma$ has the form $\widehat\gamma=5-\bar b$, where $\bar b$ estimates a posterior mean. With four hypothetical independent posterior draws $1,2,2,3$, the mean is $2$. The squared deviations sum to $2$, so $$ s_b=\sqrt{\frac{2}{3}}\approx0.8165, \qquad \operatorname{MCSE}(\bar b)=\frac{0.8165}{\sqrt{4}}\approx0.4082. $$ The update is $\widehat\gamma=5-2=3$, and its estimated simulation error is also about $0.4082$ parameter units, since subtracting $\bar b$ from the fixed value 5 does not change its standard error.

Now suppose a second batch of 16 independent draws happens to contain four copies each of $1,2,2,3$. Its mean is still $2$; its squared deviations sum to $8$, so $$ s_b=\sqrt{\frac{8}{15}}\approx0.7303, \qquad \operatorname{MCSE}(\bar b)=\frac{0.7303}{\sqrt{16}}\approx0.1826. $$ The update remains $3$, while the estimated MCSE is smaller. These hypothetical batches illustrate the calculation, not a guarantee that every larger realized batch will have a smaller estimated MCSE.

A practical schedule often increases the number of draws as the algorithm approaches convergence. Early in the iterations, changes in the parameters may be large relative to simulation noise. Later, the simulation error should be small enough that it does not obscure whether the estimated auxiliary function is improving or whether parameter updates have stabilized. There is no universal batch size that works for every model and dataset.

Useful diagnostics include checking the MCSE of important posterior expectations and of estimated changes in $Q$, repeating an E-step with a larger simulation batch, and comparing parameter updates across independent simulation batches. When comparing two values of $Q$, using the same draws to evaluate both can reduce noise in their difference. If draws are dependent, use a dependence-aware MCSE rather than treating the number of saved draws as the number of independent draws. A lack of apparent parameter movement is not persuasive evidence of convergence if the estimated update is smaller than its simulation uncertainty.

Generate posterior draws.
For each subject, simulate random effects from the posterior conditional on that subject’s observed longitudinal and time-to-event data and the current parameter value.

Build the estimated auxiliary function.
Average the complete-data log-likelihood contributions over the draws, holding the draws fixed while optimizing the estimated function.

Assess simulation uncertainty.
Estimate MCSEs for key posterior expectations or $Q$ differences, and check whether they are small relative to the changes being interpreted.

Update and reassess.
Perform the M-step, generate draws for the next iteration under the updated parameters, and increase simulation effort if Monte Carlo noise obscures progress.

Comparison, Limitations, and a Common Pitfall

Exact EM uses the posterior expectation in $Q$ without Monte Carlo approximation; when that expectation can be computed reliably, it avoids simulation error. Numerical integration within the E-step is another way to approximate the expectation and can be effective when the integration dimension and required accuracy make it feasible. MCEM is useful when the posterior expectations are too difficult to calculate accurately by those means, but its accuracy depends on the number and quality of the draws. It does not remove the need to check convergence or the fitted model.

MCEM also has a computational trade-off: improving precision requires more simulation, and the slow $M^{-1/2}$ rate can make very small MCSEs costly. An inaccurate simulation of the posterior, or draws with excessive dependence treated as independent, can make an apparently precise expectation misleading. Simulation error is additional to uncertainty arising from finite observed data; assessing one does not replace assessing the other.

Common pitfall: “Exact EM does not decrease the observed-data likelihood, so a likelihood decrease proves that my joint model or maximization is wrong.” The reasoning omits the approximation in MCEM. The exact-EM result assumes the exact posterior expectation and a valid M-step, whereas a finite Monte Carlo batch can perturb the estimated auxiliary function and its maximizer. Check the E-step simulation error and the M-step optimization before diagnosing a model or implementation error; do not assume that every decrease is harmless, either.

Key takeaway: MCEM substitutes simulated posterior averages for exact E-step expectations. Its updates can be useful and stable when the posterior draws target the right distribution and simulation error is monitored, but finite simulation does not inherit exact EM’s per-iteration likelihood guarantee.

Check Your Understanding

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

In MCEM, from which distribution must the random-effects draws for subject $i$ be generated?
If the estimated standard deviation of an independent-draw posterior quantity is $1.2$ and $M=36$, what is the estimated MCSE of its sample mean?
Why can an MCEM iteration show a decrease in observed-data likelihood even when exact EM has a non-decrease property?
A parameter update is based on a posterior expectation with MCSE $0.3$, and successive estimated updates differ by $0.05$. What diagnostic concern arises?

Solutions. (1) Draw from $p(b_i\mid D_i;\theta^{(k)})$, the subject-specific random-effects posterior conditional on observed data and current parameters. (2) $\operatorname{MCSE}=1.2/\sqrt{36}=1.2/6=0.20$ in the quantity’s units. (3) A finite simulation sample approximates the exact E-step expectation, so its estimated auxiliary function and M-step maximizer can differ from the exact ones; the exact-EM guarantee does not ensure monotonicity for each noisy MCEM iteration. (4) The estimated update is much smaller than the MCSE of the expectation on which it depends, so apparent stabilization may reflect simulation noise. Increase simulation effort or otherwise assess the uncertainty before treating the small change as evidence of convergence.