Why SAEM Needs Two Updates
In a nonlinear mixed-effects model, individual random effects are latent: they are not observed directly, but are inferred from each subject’s data and the current parameter values. SAEM, or stochastic approximation expectation-maximization, handles this by alternating between generating latent-variable values and updating estimates of complete-data sufficient statistics. It then uses those statistics in an M-step. The step size controls how quickly new simulated information replaces the current statistic estimate.
The key distinction is between the conditional distribution used to generate latent variables and the objective maximized in the M-step. A latent-variable draw is generated using the current parameter estimate. The resulting sufficient statistic is incorporated into a running estimate. The M-step then holds that running estimate fixed while optimizing over candidate parameters.
The conditional simulation may be an exact draw from \(p(z\mid y,\theta_{k-1})\), or it may use a Markov chain whose transition kernel has that conditional distribution as its invariant distribution. With Markov chain simulation, convergence arguments require suitable ergodicity and regularity conditions; a chain that does not explore the relevant conditional distribution can give persistently misleading statistics. In either case, the statistic used in the M-step is the updated \(s_k\), not a conditional expectation recomputed at every candidate value of \(\theta\).
What the Step Size Does
The recursion can also be written $$ s_k=(1-\gamma_k)s_{k-1}+\gamma_k T_k. $$ Thus \(\gamma_k\) is the weight assigned to the newest simulated statistic, while \(1-\gamma_k\) is the weight retained from the existing estimate. A large step responds quickly to new information but is more affected by the latest Monte Carlo draw. A small step smooths that noise but adapts more slowly if the conditional distribution is changing as the parameter estimates change.
A standard decreasing-step schedule is \(\gamma_k=k^{-a}\), with \(1/2<a\leq 1\). More generally, a usual stochastic-approximation condition is $$ \sum_{k=1}^{\infty}\gamma_k=\infty, \qquad \sum_{k=1}^{\infty}\gamma_k^2<\infty. $$ The first condition prevents the total amount of adaptation from being finite; the second limits the accumulated effect of random fluctuations under standard noise assumptions. These conditions are not, by themselves, a guarantee that every SAEM implementation converges. Convergence also depends on model regularity, appropriate conditional simulation, and stability of the parameter iterates.
Worked Example: Updating Random-Effect Statistics
Worked Example: Updating Random-Effect Statistics
Consider a hypothetical model with two subjects and one observation per subject: \(y_i=\theta+\eta_i+\epsilon_i\), where \(\eta_i\sim N(0,\omega^2)\), the residual errors \(\epsilon_i\) are independent with known variance \(\sigma^2\), and \(\eta_i\) and \(\epsilon_i\) are independent. The hypothetical observations are \(y_1=1.0\) and \(y_2=3.0\), in concentration units. For this model, the sufficient statistics needed to update \(\theta\) and \(\omega^2\) are \(S_1=\sum_i\eta_i\) and \(S_2=\sum_i\eta_i^2\). The residual-error variance is treated as known.
Suppose the current statistic estimates are \(s_{1,k-1}=0.40\) concentration units and \(s_{2,k-1}=0.20\) concentration-squared units. A hypothetical conditional simulation at the current parameter estimates gives \(\eta_1^{(k)}=0.20\) and \(\eta_2^{(k)}=0.80\). Therefore, $$ T_{1,k}=0.20+0.80=1.00, \qquad T_{2,k}=0.20^2+0.80^2=0.68. $$ With \(\gamma_k=0.25\), the stochastic-approximation updates are $$ s_{1,k}=0.40+0.25(1.00-0.40)=0.55, \qquad s_{2,k}=0.20+0.25(0.68-0.20)=0.32. $$
For fixed \(s_{1,k}\) and \(s_{2,k}\), the part of the M-step objective that depends on \(\theta\) and \(\omega^2\) is, up to terms constant in those candidate parameters, $$ Q_k(\theta,\omega^2) =-\frac{2}{2}\log(\omega^2)-\frac{s_{2,k}}{2\omega^2} -\frac{1}{2\sigma^2}\left[2\theta^2-2\theta\left(\sum_i y_i-s_{1,k}\right)\right]. $$ Differentiating with respect to \(\theta\) and setting the derivative to zero gives $$ \theta_k=\frac{\sum_i y_i-s_{1,k}}{2} =\frac{4.00-0.55}{2}=1.725 $$ concentration units. Maximizing with respect to \(\omega^2\) gives $$ \omega_k^2=\frac{s_{2,k}}{2}=\frac{0.32}{2}=0.160 $$ concentration-squared units. These are M-step estimates based on the updated statistics; they are not obtained by changing the conditional distribution during the maximization.
Why the Recursion Balances Adaptation and Noise
The balance can be seen first when the target statistic is fixed and the simulated statistics \(T_k\) are independent draws with mean \(m\) and variance \(v\). Subtracting \(m\) from both sides of the recursion gives $$ s_k-m=(1-\gamma_k)(s_{k-1}-m)+\gamma_k(T_k-m). $$ The first term carries forward the previous estimation error, while the second introduces fresh, mean-zero simulation noise. If \(\gamma_k\) stays large, fresh noise continues to have a substantial influence. If the steps become too small too quickly, earlier errors retain too much influence.
For a fixed target, repeated substitution expresses \(s_k\) as a weighted combination of its initial value and past draws. The weight on \(s_0\) is \(\prod_{j=1}^k(1-\gamma_j)\). When \(\sum_j\gamma_j=\infty\) and \(0<\gamma_j\leq1\), this product tends to zero: the initial estimate is eventually forgotten. Under independent, finite-variance noise, the accumulated variance from new draws is controlled by square-summable step sizes, which motivates \(\sum_j\gamma_j^2<\infty\). In SAEM, however, the target conditional expectation generally moves as \(\theta_{k-1}\) changes, and Markov chain draws can be dependent. The fixed-target argument explains the trade-off; it does not replace the additional conditions needed for convergence of the full algorithm.
Worked Example: Comparing Step Sizes During Adaptation
Suppose a hypothetical scalar sufficient statistic currently has estimate \(s_{k-1}=10.0\), and the latest simulated statistic is \(T_k=14.0\). If \(\gamma_k=0.25\), then $$ s_k=10.0+0.25(14.0-10.0)=11.0. $$ The estimate moves one quarter of the way toward the new draw. If instead \(\gamma_k=0.05\), then $$ s_k=10.0+0.05(14.0-10.0)=10.2. $$ The smaller step is less responsive to this particular draw. That can reduce the impact of a noisy draw, but it also means that a real change in the conditional target will be reflected more slowly.
Worked Example: Checking Candidate Schedules
Consider three candidate schedules, each starting at \(k=1\). For \(\gamma_k=k^{-a}\), the series \(\sum k^{-a}\) diverges when \(a\leq1\), while \(\sum k^{-2a}\) converges when \(2a>1\).
For \(\gamma_k=k^{-0.7}\), the sum of the steps diverges and the sum of their squares converges, so it meets both standard conditions. For \(\gamma_k=k^{-0.4}\), the sum of the steps diverges, but the squared-step sum \(\sum k^{-0.8}\) also diverges; this schedule does not meet the square-summability condition. For \(\gamma_k=k^{-1.2}\), the sum of the steps converges, so total adaptation is finite; it does not meet the first condition. These checks identify schedules that meet the stated conditions, not schedules guaranteed to converge regardless of model or simulation behavior.
Diagnostics, Alternatives, and Limitations
The update equations make several failure modes visible. If \(s_k\) or the parameter estimates jump sharply from iteration to iteration, inspect the size of \(\gamma_k\), the scale of the simulated statistics, and whether the conditional simulator is exploring its target distribution. If statistics change very little despite evidence that the target is shifting, the schedule may be shrinking too quickly. A sufficient-statistic estimate that becomes invalid for the intended parameterization, such as a negative estimated variance when the M-step requires a positive variance, signals a problem with the update, the chosen statistic, or the implementation. The M-step must respect the model’s parameter space; an unconstrained formula cannot be used blindly if its result is inadmissible.
SAEM differs from the likelihood approximations covered in the earlier material on FO, FOCE, and Laplace: it uses conditional simulation and stochastic approximation of complete-data statistics rather than relying on those approximations. Those alternatives may be preferable when their assumptions and computational behavior suit the model and design. SAEM can be useful when latent-effect integrals are difficult to evaluate directly, but it introduces Monte Carlo variability and depends on a functioning conditional simulator. It also requires an M-step that can be carried out for the chosen statistics; the exponential-family representation makes that step explicit, but not every model yields a simple closed-form maximizer.
A common mistake is to define the M-step objective using a conditional distribution that changes with each candidate parameter value. That would make the latent-data expectation vary during maximization. In SAEM, the conditional distribution for generating the iteration’s latent data uses the pre-M-step estimate \(\theta_{k-1}\); the resulting statistic estimate \(s_k\) is then held fixed while maximizing \(Q_k(\theta)\). This distinction is essential to the update sequence.
Using the Updates in a Reproducible Workflow
Identify the latent variables and sufficient statistics, and state the parameter constraints and assumptions used by the model.
At iteration \(k\), use the conditional distribution associated with \(y\) and \(\theta_{k-1}\), through a valid draw or a suitable Markov chain transition.
Form \(T_k\) from the simulated latent data and apply \(s_k=s_{k-1}+\gamma_k(T_k-s_{k-1})\).
Compute \(\theta_k\) from \(Q_k(\theta)=\langle\phi(\theta),s_k\rangle-A(\theta)\), respecting parameter constraints.
Check that the schedule is appropriate, the statistics remain plausible, and changes in estimates are consistent with the amount of adaptation allowed.
For reproducibility, record the step-size formula and its constants, the initialization of \(s_0\) and \(\theta_0\), and the random-number seed and simulation settings used to generate latent values. These details make it possible to distinguish a change in model behavior from a change in the stochastic update process.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions below.
- At iteration \(k\), which parameter values define the conditional distribution used to generate latent data, and which statistic is held fixed during the M-step?
- If \(s_{k-1}=6\), \(T_k=10\), and \(\gamma_k=0.30\), calculate \(s_k\).
- Does \(\gamma_k=k^{-0.6}\) satisfy the two standard step-size sum conditions? Explain.
- Why can a fixed step size leave ongoing Monte Carlo variability in the statistic estimate?
- Does satisfying the two step-size sum conditions alone guarantee convergence of a particular SAEM fit? Why or why not?
Solutions. (1) The latent data are generated conditional on \(y\) and \(\theta_{k-1}\); the updated statistic \(s_k\) is held fixed while candidate \(\theta\) values are considered in the M-step. (2) \(s_k=6+0.30(10-6)=7.2\). (3) Yes. Since \(0.6\leq1\), \(\sum k^{-0.6}\) diverges; since \(1.2>1\), \(\sum k^{-1.2}\) converges. (4) A fixed positive step continues to give new simulated values nonzero weight, so their random variation continues to affect the update. (5) No. Appropriate conditional simulation, model regularity, and stability of the parameter iterates are also needed.