Why SAEM Uses a Running Approximation
In a nonlinear mixed-effects model, each subject has an unobserved individual effect. The complete-data likelihood is often easier to work with than the likelihood obtained after integrating those effects out, but its EM update requires a conditional expectation over the effects. In nonlinear models that expectation may not have a convenient exact calculation. SAEM replaces the exact conditional expectation with a simulation-based, recursively updated approximation.
The key is that SAEM does not simply maximize the complete-data likelihood for one simulated set of effects. A single draw can be noisy. Instead, each iteration blends newly simulated complete-data statistics with the running estimate from earlier iterations, then maximizes the resulting approximate criterion.
The SAEM Iteration
The standard presentation assumes the complete-data log-likelihood can be written, up to a term not depending on \(\theta\), in terms of a statistic vector \(S(\mathbf{y},\boldsymbol{\eta})\):
Here \(c\) does not depend on \(\theta\), \(\psi\) is a scalar function, \(\phi(\theta)\) is a parameter-dependent vector, and \(\langle\cdot,\cdot\rangle\) is the vector inner product. The particular form of \(S\) depends on the model; the essential point for this tutorial is that the complete-data criterion depends on the latent effects through this statistic.
At iteration \(k\), SAEM has a current parameter value \(\theta_k\) and a running statistic \(s_{k-1}\). It performs three operations:
Generate \(\boldsymbol{\eta}^{(k)}\) from the conditional distribution of latent effects given the observed data at the current parameter value, \(p(\boldsymbol{\eta}\mid\mathbf{y},\theta_k)\). In practice, the simulation method must target this conditional distribution adequately.
Calculate \(S_k^{\mathrm{sim}}=S(\mathbf{y},\boldsymbol{\eta}^{(k)})\) and blend it with the previous statistic using a gain \(\gamma_k\).
Use the updated statistic \(s_k\) in place of the unavailable exact conditional expectation, and maximize the resulting criterion to obtain \(\theta_{k+1}\).
The update in Step 2 is \(s_k=s_{k-1}+\gamma_k(S_k^{\mathrm{sim}}-s_{k-1})\), with \(0<\gamma_k\leq 1\). The corresponding approximate criterion is \(Q_k(\theta)=-\psi(\theta)+\langle s_k,\phi(\theta)\rangle\), up to terms that do not depend on \(\theta\). Thus the parameter update is \(\theta_{k+1}\in\arg\max_\theta Q_k(\theta)\), when a maximizer exists.
The draws at iteration \(k\) use \(\theta_k\), and the maximization produces \(\theta_{k+1}\). Keeping that indexing straight matters: the parameter update is based on the statistic updated using draws from the current conditional distribution, not from a distribution that already uses \(\theta_{k+1}\).
Why the Recursion Approximates the EM Step
For fixed \(\theta_k\), the conditional mean of the simulated statistic is \( \bar S(\theta_k)=E[S(\mathbf{y},\boldsymbol{\eta})\mid\mathbf{y},\theta_k] \). The exact EM criterion, apart from terms not depending on the candidate parameter, is \( -\psi(\theta)+\langle\bar S(\theta_k),\phi(\theta)\rangle \). SAEM substitutes its running estimate \(s_k\) for \(\bar S(\theta_k)\).
The reason a running estimate can be useful follows directly from the recursion. Conditional on the current parameter and the past, suppose the simulated statistic has mean \(\bar S(\theta_k)\). Then:
So the expected update moves a fraction \(\gamma_k\) of the gap from the current running estimate toward the conditional mean. A small gain dampens the influence of an individual draw; a gain of one replaces the running estimate with the newest simulated statistic. This is the stochastic-approximation mechanism, not an exact evaluation of the E-step.
For a changing parameter value, the conditional mean also changes over iterations. Standard convergence arguments therefore require conditions on the model and the simulation as well as a gain sequence. A common set of gain conditions is \(0<\gamma_k\leq1\), \(\sum_k\gamma_k=\infty\), and \(\sum_k\gamma_k^2<\infty\). These conditions describe a sequence that continues to respond to the target while reducing the cumulative effect of simulation noise. They do not by themselves guarantee convergence if the conditional simulation is inaccurate, the likelihood is irregular, or the M-step is not carried out reliably.
Worked Updates
Worked Example: One Simulated Statistic at Each Iteration
Consider a hypothetical complete-data criterion with scalar statistic \(T\): \( \ell_c(\theta)=\theta T-\theta^2/2+c \), where \(c\) does not depend on \(\theta\). Its maximizer is \(\theta=T\), since differentiating gives \(T-\theta=0\), and the second derivative is \(-1\). Suppose simulated values of \(T\) at three successive iterations are \(2.0\), \(0.0\), and \(1.5\). Start with \(s_0=0\), and use gains \(\gamma_1=0.5\), \(\gamma_2=0.25\), and \(\gamma_3=0.2\).
At the first iteration, \(s_1=0+0.5(2.0-0)=1.0\), so the maximizer of the approximate criterion is \(\theta_2=1.0\).
At the second iteration, \(s_2=1.0+0.25(0.0-1.0)=0.75\), giving \(\theta_3=0.75\).
At the third iteration, \(s_3=0.75+0.2(1.5-0.75)=0.75+0.15=0.90\), giving \(\theta_4=0.90\).
Each parameter update maximizes the criterion based on the current running statistic. The sequence of simulated values itself is not the sequence of parameter estimates: the gains determine how much each new value changes the statistic used in the M-step.
Worked Example: Using Conditional Latent-Effect Draws
Suppose that, at one SAEM iteration, the simulated latent effects for a hypothetical subject yield scalar statistics \(0.8\), \(-0.4\), and \(0.5\). Treat these as three draws from the conditional distribution at the current parameter value, and average them for this iteration: \(S_k^{\mathrm{sim}}=(0.8-0.4+0.5)/3=0.9/3=0.3\). Assume the running statistic before this update is \(s_{k-1}=0.1\), and the gain is \(\gamma_k=0.4\).
The stochastic-approximation update is \(s_k=0.1+0.4(0.3-0.1)=0.1+0.08=0.18\). If the approximate criterion for the population parameter is the one in the previous example, its maximizer is \(\theta_{k+1}=s_k=0.18\).
The sample average \(0.3\) is the statistic used in the recursion for this iteration; it is not asserted to equal the exact conditional expectation. Repeating the procedure can reduce the effect of simulation variability through the running update. The calculation presumes the draws are generated from the intended conditional distribution; draws from an unsuitable distribution would update the statistic toward the wrong target.
Worked Example: Combining Statistics Across Subjects
For a hypothetical model with three subjects, suppose the complete-data criterion for a scalar parameter is \( \ell_c(\theta)=\theta(T_1+T_2+T_3)-3\theta^2/2+c \). The sufficient statistic is the sum \(T_1+T_2+T_3\), and maximizing this criterion gives \(\theta=(T_1+T_2+T_3)/3\). At an iteration, simulated subject-level statistics are \(0.2\), \(0.8\), and \(-0.1\), so their sum is \(0.9\). Let the previous running sum be \(0.3\) and the gain be \(0.5\).
Update the sum: \(s_k=0.3+0.5(0.9-0.3)=0.3+0.3=0.6\). Because the criterion uses the sum and there are three subjects, the M-step is \(\theta_{k+1}=s_k/3=0.6/3=0.20\).
The division by three is part of this particular criterion’s maximization. A statistic should not be divided or averaged by habit: the M-step must follow from the actual parameter-dependent complete-data criterion.
Assumptions, Diagnostics, and Alternatives
SAEM’s interpretation relies on several conditions. The model must define a valid joint distribution for the observations and latent effects, and the conditional distribution used for simulation must correspond to \(p(\boldsymbol{\eta}\mid\mathbf{y},\theta_k)\). The chosen statistic must retain the dependence of the complete-data criterion on the population parameters. The M-step must maximize the approximate criterion, or at least increase it sufficiently for the chosen implementation. For convergence arguments, the statistics should have suitable finite moments, the gain sequence should meet appropriate stochastic-approximation conditions, and the parameterization should not create uncontrolled boundary or identifiability problems.
These conditions suggest practical checks. Monitor the running statistics and parameter estimates across iterations: persistent large jumps or drift can indicate excessive simulation noise, a gain that does not diminish appropriately, or an unstable maximization. Check that the conditional simulation explores the relevant latent-effect regions rather than repeatedly returning nearly identical values for computational reasons. Re-evaluate the approximate criterion after the M-step and inspect whether optimization has actually converged. Compare fits from different initial values, since convergence to a stationary point does not establish that it is the global optimum.
A particularly common mistake is to treat one simulated latent-effect sample as if it were the exact E-step. With \(\gamma_k=1\), the recursion does exactly replace the running statistic with that iteration’s simulated statistic; that is not generally an accurate conditional expectation. The correction is to recognize the role of the gain and, when appropriate, use a sequence that reduces the influence of noisy draws while still allowing the estimate to track its changing target.
Exact EM instead computes the conditional expectation directly, which can be difficult in nonlinear mixed models. Monte Carlo EM estimates that expectation using simulation for an iteration, whereas SAEM recursively smooths simulated statistics across iterations. Direct maximization of the observed-data likelihood avoids this EM-style latent-statistic update but requires its own reliable likelihood calculations and optimization. These approaches have different computational demands; none removes the need to assess numerical stability and the adequacy of the fitted model.
SAEM can be sensitive to poor conditional simulation, weakly identified parameters, and an ill-conditioned M-step. Its stochastic updates can obscure progress in short runs, while an overly persistent gain can leave substantial Monte Carlo variability. Conversely, decreasing the gain too aggressively before the statistics have approached their relevant targets can make later updates ineffective. The method is an estimation algorithm, not a diagnostic for whether the model’s structural or residual-error assumptions are appropriate.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- If the simulated statistic is \(1.4\), the previous running statistic is \(0.6\), and \(\gamma_k=0.25\), what is the updated statistic?
- Why does SAEM simulate latent effects conditional on both the observed data and the current parameter value?
- In a complete-data criterion \(\ell_c(\theta)=\theta T-\theta^2/2+c\), what parameter maximizes the criterion when the running statistic is \(s_k=0.35\)?
- What does setting \(\gamma_k=1\) do to the running statistic at that iteration?
- Name one reason that apparently stable parameter estimates would not, by themselves, establish that SAEM has converged to the global maximum.
Solutions. (1) \(s_k=0.6+0.25(1.4-0.6)=0.6+0.2=0.8\). (2) That conditional distribution is the distribution of the latent effects given the data under the current model parameters; it supplies the conditional complete-data statistics needed for the approximate E-step. (3) The approximate criterion is \(\theta(0.35)-\theta^2/2+c\), whose derivative is \(0.35-\theta\), so its maximizer is \(\theta=0.35\). (4) The update becomes \(s_k=s_{k-1}+(S_k^{\mathrm{sim}}-s_{k-1})=S_k^{\mathrm{sim}}\): it replaces the running statistic with the new simulated statistic. (5) SAEM may converge to a local stationary point, so different initial values or other evidence may be needed to assess whether a better maximum exists.