A Worked SAEM Fit for Nonlinear Exponential Decay

Why Work Through a Nonlinear SAEM Fit?

In an exponential-decay mixed model, each subject has an individual amplitude and decay rate, while the population distribution describes how those subject-specific quantities vary. The conditional mean is nonlinear in the effects, so the conditional expectations required by exact EM are generally not available in closed form. SAEM replaces those expectations with simulated latent effects and smooths the resulting sufficient statistics across iterations.

The goal here is to follow one numerical iteration from simulated effects to updated population parameters. The example is deliberately small and hypothetical; it illustrates the arithmetic, not a real pharmacometric analysis. We assume the latent effects have been generated from their conditional distribution given the observations and current parameters. The method used to generate those draws is a separate implementation issue.

Model and Complete-Data Quantities

There are \(m\) independent subjects. Subject \(i\) is observed at times \(t_{ij}\), for \(j=1,\ldots,n_i\), with total observation count \(N=\sum_i n_i\). Write \(a_i\) for log amplitude and \(b_i\) for log decay rate. The model is

The random-effects distribution is independent across subjects, and within a subject we assume the two effects are independent:

The population parameter vector is \(\theta=(\mu_a,\mu_b,\omega_a^2,\omega_b^2,\sigma^2)\). Conditional on the effects, observations are independent normal variables with the stated means and common residual variance. Effects are independent across subjects, and the random-effects covariance is diagonal. These are modeling assumptions, not consequences of SAEM.

Ignoring constants that do not depend on the parameters, the complete-data log-likelihood is

The nonlinear observation model enters through \(Q(a,b)\). SAEM approximates the conditional expectations of the complete-data quantities using simulated effects, then updates smoothed statistics rather than replacing them with each iteration’s raw simulation.

Stochastic Approximation and the M-Step

Let \(s^{(r)}\) be the vector of complete-data quantities calculated from the latent-effect draw at iteration \(r\). Let \(S^{(r)}\) be their running stochastic-approximation values. With gain \(\gamma_r\), the update is

When \(\gamma_r=1\), the statistic is replaced by the latest simulated value. With a gain between zero and one, the update moves only partway toward that value. The M-step maximizes the expected complete-data log-likelihood represented by the smoothed statistics.

For the Gaussian random effects, write the running statistics as \(S_a,S_{a2},S_b,S_{b2}\). Their M-step updates are

If \(S_Q\) is the running statistic for the residual sum of squares, then \(\widehat{\sigma}^2=S_Q/N\). These formulas presume the diagonal covariance and common residual variance specified above. In an actual implementation, a variance estimate at or below zero signals a numerical or modeling problem; it is not a valid Gaussian variance.

One SAEM Iteration, Calculated

Why Update Statistics Rather Than Parameters Directly?

The stochastic approximation applies to the complete-data statistics. The parameter update is then computed from those smoothed statistics. Smoothing parameter estimates themselves is generally not the same operation: the M-step is nonlinear, as the variance formula’s subtraction of a squared mean makes clear.

The reason the formulas have this form follows directly from the complete-data likelihood. For \(a_i\), the part involving the population mean and variance is proportional to \(-m\log(\omega_a^2)/2-\sum_i(a_i-\mu_a)^2/(2\omega_a^2)\). Differentiating with respect to \(\mu_a\) and setting the result to zero gives \(\sum_i(a_i-\mu_a)=0\), hence \(\mu_a=\sum_i a_i/m\). Substituting that mean and maximizing with respect to \(\omega_a^2\) gives \(\omega_a^2=\sum_i(a_i-\mu_a)^2/m\), which expands to \(\sum_i a_i^2/m-\mu_a^2\). Replacing the complete-data sums by their smoothed SAEM statistics yields the displayed updates. The same argument applies to \(b_i\). Differentiating the residual part \(-N\log(\sigma^2)/2-S_Q/(2\sigma^2)\) gives \(\sigma^2=S_Q/N\).

Diagnostics, Alternatives, and Limitations

A single SAEM iteration is not expected to improve every statistic monotonically. Simulation noise can make running estimates move up or down, and the observed-data likelihood need not increase at every stochastic iteration. Assess the fit over a run, not from one update alone.

• Check the simulated effects. Draws that barely change across iterations may indicate inadequate exploration of the conditional distribution. A narrow-looking trace is not reassuring if the sampler is stuck.
• Inspect the statistics and parameter traces. Persistent large oscillations can indicate excessive simulation noise, an unsuitable gain, or weak information in the data. A variance approaching zero or an implausibly large value warrants investigation.
• Check the model fit over time. Residual patterns across observation time can reveal a poor exponential-decay description, even when the residual variance update is numerically valid.
• Assess information in the design. Sparse observations, especially observations over a narrow time range, can make amplitude and decay rate difficult to distinguish. SAEM cannot create information that the sampling design does not contain.

Exact EM uses conditional expectations rather than simulated approximations, but those expectations are often intractable for nonlinear random-effects models. Monte Carlo EM estimates the expectations afresh by simulation at each iteration; SAEM instead recursively smooths sufficient statistics, which can reduce the need for highly accurate Monte Carlo averages at every step. Direct likelihood maximization with numerical integration is another approach, but its reliability depends on accurate integration. These alternatives have different computational and numerical trade-offs; none removes the need to check model fit and identifiability.

The example assumes independent random effects and constant additive residual variance. It does not address correlated effects, alternative residual-error models, or the details of constructing a conditional simulator. Those choices can change the sufficient statistics, M-step, or simulation behavior. In particular, the estimates above are one iteration’s updates and should not be reported as a converged fit.