Common SAEM Implementation Mistakes

Why SAEM Implementations Need Auditing

A SAEM iteration combines conditional simulation, a stochastic-approximation update, and parameter maximization. Each step can look numerically plausible while using the wrong quantity: a gain may be applied to the new statistic rather than to its difference from the current estimate; simulated latent effects may correspond to an unintended parameter vector; or a covariance matrix may cease to be a valid covariance. These errors can produce smooth-looking traces without implementing the intended algorithm.

The aim here is not to restate the general theory of SAEM. It is to make the iteration auditable: record what parameter vector generated each conditional simulation, verify the exact gain formula, and check the covariance produced by each update. The examples use hypothetical values solely to illustrate calculations.

The Iteration to Be Checked

Let \(y\) denote the observed data, \(z\) the collection of latent effects or other latent variables, and \(\theta\) the model parameter vector. Write \(s(y,z)\) for the complete-data sufficient-statistic vector used by the M-step. In a typical SAEM iteration, a latent draw \(z_k\) is simulated from its conditional distribution given the observations and the current parameter vector, a running statistic is updated, and the parameters are maximized using that statistic.

The statistic update is an exponential moving average: it retains weight \(1-\gamma_k\) on the old estimate and assigns weight \(\gamma_k\) to the new simulated value. In standard stochastic-approximation schedules, gains are often chosen to decrease while satisfying \(\sum_k\gamma_k=\infty\) and \(\sum_k\gamma_k^2<\infty\), subject to the algorithm’s regularity conditions. Implementations may also use an initial finite phase with gain one. Whatever schedule is chosen, the code must use the gain attached to the intended iteration and apply it consistently.

Mistake 1: Applying the Gain to the Wrong Quantity

The most direct implementation error is to code \(S_k=S_{k-1}+\gamma_k s(y,z_k)\). This adds a fraction of the new statistic to the old statistic but never removes the corresponding fraction of the old value. It is not the stated stochastic-approximation update. The correct algebraic form is either the definition above or its equivalent weighted-average form:

Other gain mistakes include using \(\gamma_{k-1}\) instead of \(\gamma_k\), applying a gain twice (once in a helper function and again in the calling code), or resetting \(S_{k-1}\) at every iteration. To detect them, log the gain, old statistic, simulated statistic, and updated statistic together, and independently recompute the update for selected iterations. For \(0\leq\gamma_k\leq1\), a scalar update must lie between \(S_{k-1}\) and \(s(y,z_k)\); a value outside that interval signals an error.

Mistake 2: Simulating Under an Unintended Parameter Snapshot

A conditional simulation is defined relative to a full parameter vector. If the algorithm prescribes simulation using \(\theta_{k-1}\), then the likelihood and latent-effect distribution used by the conditional sampler must both be evaluated at that snapshot. A component-wise mixture of old and new values can define a perfectly valid conditional distribution; the issue is whether it is the distribution required by the stated iteration.

In a nonlinear mixed-effects model, the conditional distribution may be sampled by an MCMC transition or another numerical procedure rather than obtained from a normal formula. The same audit applies: identify the parameter snapshot defining the target, ensure every model component used by the conditional calculation corresponds to that snapshot, and do not silently replace components partway through a simulation sweep. Any intentional alternative scheduling must be stated as part of the algorithm and implemented consistently.

A parameter change after simulation is not itself an error: in the iteration defined above, the M-step produces \(\theta_k\) after the conditional draw. The key distinction is that the simulation and M-step occupy different stages and therefore can use different parameter values. A trace that stores only the final \(\theta_k\) may conceal this distinction, so record both the simulation snapshot and the updated parameter vector.

Mistake 3: Producing an Invalid Covariance Update

A covariance matrix must be symmetric and positive semidefinite; for a nondegenerate multivariate normal random-effect distribution, it must be positive definite. For a two-dimensional covariance matrix \(\Omega=\begin{pmatrix}v_1&c\\c&v_2\end{pmatrix}\), these conditions include \(v_1\geq0\), \(v_2\geq0\), and \(v_1v_2-c^2\geq0\). Positive definiteness requires strict inequalities for the eigenvalues, equivalently positive diagonal entries and positive determinant in the two-dimensional case.

When the random effects have mean zero and covariance \(\Omega\), the covariance contribution to the Gaussian complete-data log likelihood depends on the second moment \(bb^\mathsf{T}\). Consequently, a typical covariance M-step uses the average of the corresponding expected second moments, not a sample covariance centered by an estimated sample mean unless that is explicitly part of the model’s parameterization. This distinction matters when translating sufficient statistics into a covariance update.

There is also a useful justification for checking a matrix-valued stochastic-approximation update. If \(A\) and \(B\) are positive semidefinite and \(0\leq\gamma\leq1\), then \(C=(1-\gamma)A+\gamma B\) is positive semidefinite because for every vector \(x\),

Thus a convex update of valid second-moment matrices preserves positive semidefiniteness in exact arithmetic. This does not guarantee positive definiteness: a convex combination can remain singular, for example if both inputs share a null direction. Nor does it rescue an invalid input statistic, an incorrectly assembled asymmetric matrix, or an M-step that does not correspond to the stored statistics.

A Practical Audit and Its Limits

These checks complement, rather than replace, convergence diagnostics. A stable trace cannot prove that the intended algorithm was implemented: a faulty gain or unintended conditional target may also yield a stable sequence. Conversely, occasional noisy movement is not by itself evidence of a coding defect in a stochastic algorithm. Diagnose the algebra and parameter flow first, then assess Monte Carlo variability and convergence under the chosen method.

The main alternatives do not remove the need for care. Monte Carlo EM typically estimates an expectation using multiple conditional draws, whereas SAEM updates running statistics with a gain; confusing those update mechanisms changes the algorithm. Deterministic EM avoids simulation noise when its conditional expectations are tractable, but it does not apply to every nonlinear model. In all cases, the target distribution, statistic definitions, and covariance parameterization must match the model. SAEM’s flexibility also has limits: finite simulation, poor mixing, weak identification, or an unsuitable parameterization can impair estimation even when the code follows the update equations exactly.

Common Mistakes to Avoid

• Adding the new statistic without subtracting the old one. The correct increment is \(\gamma_k(s_k-S_{k-1})\), not \(\gamma_k s_k\).
• Calling every hybrid conditional “invalid.” A fully specified hybrid parameter vector can define a valid distribution; it is an implementation error only when it is not the target prescribed for that stage.
• Checking variances but not the whole covariance. Positive diagonal entries do not guarantee a positive-semidefinite matrix.
• Using a successful run as proof of correctness. Convergence-looking traces do not verify the gain formula, parameter snapshot, or M-step inputs.
• Repairing a bad covariance without investigating it. A numerical adjustment may conceal a statistic or coding error; first identify why the proposed update violates the covariance constraints.