From a Model to an Estimation Cycle
SAEM combines conditional simulation of latent random effects with stochastic approximation of complete-data statistics. An end-to-end calculation is useful because these pieces must agree: the conditional simulation uses the current model, the gain controls how much a new simulated state changes the running statistics, and the M-step updates parameters using those statistics. A plausible parameter estimate is not enough if the steps that produced it are inconsistent.
Consider a hypothetical decay study with \(m=3\) subjects and two response measurements per subject, at \(t=1\) and \(t=2\) hours. Let \(y_{ij}\) denote the response for subject \(i\) at time \(t_{ij}\). The response is dimensionless. Use the nonlinear model
The random effects \(b_i\) are independent across subjects; the residuals are independent across observations and independent of the random effects. The population parameter \(\theta\) is a log-rate, so \(\exp(\theta)\) has units of inverse hours. Here \(\omega^2\) is the random-effect variance on the log-rate scale, and \(\sigma^2\) is the residual variance on the response scale. Conditional on \(b_i\), the observations are normally distributed about the stated nonlinear mean.
Complete-Data Criterion and the SAEM Update
Write \(\psi=(\theta,\omega^2,\sigma^2)\), and define the model mean \(\mu_{ij}(\theta,b_i)=\exp\{-t_{ij}\exp(\theta+b_i)\}\). Ignoring constants that do not depend on \(\psi\), the complete-data log likelihood is
where \(N=\sum_i n_i\) is the number of observations. The normal random-effects distribution has mean zero, as specified by the model; it is not estimated as an additional mean parameter in this example. For fixed \(\theta\), the complete-data criterion depends on the latent effects through \(\sum_i b_i^2\) and the residual sum of squares. Because the mean is nonlinear in \(\theta\), maximizing the criterion over \(\theta\) generally requires numerical optimization.
Let \(S_k\) denote an accumulated statistic and \(s_k\) the corresponding statistic calculated from the current conditional simulation. A standard stochastic-approximation update is
The gain \(\gamma_k\) is between zero and one. Early iterations may use larger gains to adapt quickly; later gains are commonly reduced, with conditions such as \(\sum_k\gamma_k=\infty\) and \(\sum_k\gamma_k^2<\infty\) used in standard convergence arguments. A finite implementation must also use a suitable conditional simulation, adequate exploration, and a well-behaved M-step. The gain conditions alone do not guarantee a correct result.
Worked Example: One Conditional Random-Effect Move
Worked Example: One Conditional Random-Effect Move
For one hypothetical subject, the observations are \(y(1)=0.70\) and \(y(2)=0.50\). At the current values \(\theta=-1.00\), \(\omega^2=0.0400\), and \(\sigma^2=0.0025\), consider a symmetric random-walk proposal from \(b=0\) to \(b^\star=-0.10\). The conditional target for \(b\) is proportional to the subject’s likelihood times its \(N(0,\omega^2)\) density. A symmetric proposal has cancelling forward and reverse proposal densities, so the log acceptance ratio is
At \(b=0\), the rate is \(\exp(-1)=0.3679\) per hour. The model predictions are \(\exp(-0.3679)\approx0.6922\) and \(\exp(-0.7358)\approx0.4791\). Thus
At \(b^\star=-0.10\), the rate is \(\exp(-1.10)=0.3329\) per hour, and the predictions are approximately \(0.7169\) and \(0.5139\). Therefore
Using the unrounded model values, the likelihood contribution to \(\log r\) is about \(- (0.000477334-0.000495899)/(0.005)=0.003713\). The random-effect contribution is \(-0.0100/0.0800=-0.1250\), giving \(\log r\approx-0.121\) and \(r\approx0.886\). The acceptance probability is \(\min(1,r)=0.886\). If a uniform draw is below that probability, accept \(b^\star\); otherwise retain \(b=0\). This is a move in the conditional simulation, not a claim that the proposed value is more likely under the prior alone.
Worked Example: Stochastic Approximation and the M-Step
Worked Example: Updating Statistics and Variances
Suppose the current conditional simulation gives the three subject effects \(b=(-0.10,0.18,0.06)\). Their raw squared-effect statistic is
For illustration, take the observations to be the model means at \(\theta=-1\) and these effects, plus residuals chosen as follows. For each subject, let \(e_{i1}=c_i\), where \(c=(-0.010,0.020,-0.010)\), and set \(e_{i2}=-e_{i1}w_{i1}/w_{i2}\), with \(w_{ij}=t_{ij}\exp(\theta+b_i)\mu_{ij}(\theta,b_i)\). This defines the hypothetical observations by \(y_{ij}=\mu_{ij}+e_{ij}\). It also makes the derivative of the residual sum of squares with respect to \(\theta\) equal to zero at \(-1\), since \(\partial\mu_{ij}/\partial\theta=-w_{ij}\) and \(\sum_{ij}e_{ij}w_{ij}=0\). The residual sum of squares at that value is approximately \(0.0009446\), or \(s_e=0.0009446\) for this simulated complete-data state. A one-dimensional numerical check around \(-1\) confirms a local minimum of the current state's residual sum of squares, but this alone does not establish the maximizer of the accumulated SAEM objective, which also depends on previously accumulated statistics.
Suppose the previous accumulated statistics were \(S_{b,k-1}=0.0200\) and \(S_{e,k-1}=0.00110\), and the gain is \(\gamma_k=0.25\). The updated statistics are
For \(m=3\) effects and \(N=6\) observations, maximizing the normal complete-data terms gives \(\widehat{\omega}^2=S_{b,k}/m=0.00883\) and \(\widehat{\sigma}^2=S_{e,k}/N=0.000177\). These are variance estimates on their respective scales. In a practical nonlinear M-step, the residual statistic used in the objective must be represented consistently as a function of candidate \(\theta\); one must not optimize \(\theta\) using residuals frozen at an unrelated parameter value. The current-state stationarity at \(\theta=-1\) does not determine the M-step update; that requires optimizing the accumulated objective as a function of candidate \(\theta\).
Why the M-Step and Convergence Checks Matter
The variance formulas above follow by differentiating the relevant normal terms. For example, the effect-dependent term is \(-m\log(\omega^2)/2-S_b/(2\omega^2)\); its derivative with respect to \(\omega^2\) is \(-m/(2\omega^2)+S_b/(2(\omega^2)^2)\), which is zero at \(\omega^2=S_b/m\). The residual variance follows in the same way, replacing \(m,S_b\) with \(N,S_e\). The update for \(\theta\) has no comparable closed form here: its objective contains exponentials of \(\theta+b_i\), so a numerical maximizer is required. Positive variance estimates and a positive-definite covariance are necessary; in this one-effect example, positive variance means simply \(\omega^2>0\).
In a full SAEM run, those updates repeat. For example, conditional simulations may use a Metropolis–Hastings kernel targeting \(p(b\mid y,\psi)\); after a burn-in or suitable within-iteration simulation, the resulting complete-data statistics enter the gain update. The same parameterization and likelihood must be used throughout. The stochastic update is not a direct average of parameter estimates: it averages the chosen statistics, after which the M-step maps them to parameters.
A convergence assessment should consider parameter traces and accumulated statistics together. Stable-looking parameter estimates with continuing drift in sufficient statistics may indicate that the run has not stabilized. Conversely, noisy individual simulations do not by themselves show parameter instability if decreasing gains are producing stable statistics and estimates. Inspect acceptance behavior and the range of simulated effects as well: extremely low acceptance can leave chains nearly fixed, while poor exploration can make the accumulated statistics unrepresentative of the conditional distribution. Compare independent runs from reasonable initial values when local maxima or initialization sensitivity are plausible. These checks support a numerical assessment; they do not prove that a global maximum has been found.
Worked Example: Reading a Late-Run Trace
Suppose a hypothetical run reports the following values at iterations 400, 450, and 500: \(\theta=(-1.012,-1.009,-1.011)\), \(\omega^2=(0.082,0.081,0.082)\), and \(\sigma^2=(0.0021,0.0021,0.0020)\). Over the same iterations, the accumulated residual statistic changes by less than one percent, and the conditional random-effect sampler accepts about 30% of proposals without a sustained trend in its simulated effects. The small parameter movements and stable statistics are evidence of late-run stabilization, not a formal guarantee of convergence. If the statistic continued to trend while the parameters appeared nearly constant, the run should continue or be investigated rather than declared converged on the parameter trace alone.
SAEM is attractive when the latent-effects conditional distribution cannot be integrated conveniently and a useful complete-data structure yields manageable M-steps. Monte Carlo EM instead uses a Monte Carlo approximation to the conditional expectation at each iteration; its simulation error can obscure small changes in the EM objective unless simulation effort increases. Deterministic integration can be effective for low-dimensional random effects, but may become computationally difficult as their dimension grows. SAEM still depends on an adequate conditional simulator and can be sensitive to scaling, gain schedules, initialization, and multimodality. Its stabilized estimates are maximum-likelihood estimates, not posterior summaries.
Common Mistake: Treating One Cycle as Convergence
A common mistake is to interpret one accepted random-effect proposal and one M-step as a completed estimation. That confuses a single stochastic update with a converged likelihood optimization. SAEM’s conditional simulation is noisy, and a single draw need not approximate a conditional expectation well. The correction is to iterate with a suitable gain schedule, monitor statistics as well as parameters, check simulation behavior, and assess sensitivity to initialization. Also keep the scales straight: \(b_i\) is dimensionless here, \(\omega^2\) is a variance on the log-rate scale, and \(\sigma^2\) is a variance in squared response units.
Check Your Understanding
Use the model and definitions above. The brief solutions follow the questions.
- Why does a symmetric random-walk proposal cancel from the conditional Metropolis–Hastings ratio?
- If \(S_{k-1}=0.030\), \(s_k=0.050\), and \(\gamma_k=0.20\), what is \(S_k\)?
- With three subjects and an accumulated squared-effect statistic of \(0.0265\), what is the corresponding M-step estimate of \(\omega^2\)?
- Why is a stable parameter trace alone insufficient evidence that SAEM has stabilized?
Solutions. (1) Symmetry gives equal forward and reverse proposal densities, so their ratio is one. (2) \(S_k=0.80(0.030)+0.20(0.050)=0.034\). (3) \(\widehat{\omega}^2=0.0265/3=0.00883\), on the log-rate random-effect scale. (4) The statistics may still drift, or the conditional sampler may explore poorly; inspect both statistics and parameters, along with simulation behavior.