Closed-Form M-Steps Within SAEM

Why Some SAEM M-Steps Have Closed-Form Updates

A stochastic approximation expectation-maximization (SAEM) iteration can involve difficult conditional expectations, particularly when the population model is nonlinear. But once SAEM has accumulated the sufficient statistics needed for a parameter block, maximizing the corresponding expected complete-data log likelihood may be straightforward. In particular, Gaussian variance parameters often have direct updates that require only a sum of squared latent effects or residuals and its associated count.

The key is to keep separate the two stages: SAEM estimates the required conditional sufficient statistics, then the M-step maximizes the relevant likelihood expression using those estimates. The updates below do not require replacing the nonlinear population model with a linear one. They do require the variance parameters to enter the complete-data likelihood in the specified Gaussian form, and they are conditional on any mean or nonlinear parameters being held fixed or updated in a separate, compatible step.

Let \(L=\sum_{i=1}^N n_i\) be the total number of observations. At an SAEM iteration, denote the current approximations to the relevant conditional sufficient statistics by \(S_{b^2}\), an estimate of \(\sum_i E[b_i^2\mid y]\), and \(S_{e^2}(\phi)\), an estimate of \(\sum_{i,j}E[(y_{ij}-f_{ij}(\phi,b_i))^2\mid y]\). The expectations are conditional on the observed data, with the latent-effect distribution evaluated under the current E-step parameters. In practice, SAEM updates approximations to such statistics using simulated latent effects and stochastic approximation; these \(S\) quantities are not necessarily the squared values from one simulated draw.

The notation \(S_{e^2}(\phi)\) makes an important dependency explicit: if the prediction parameters \(\phi\) change, the residuals in this statistic can change too. The residual variance update is closed form when the residual-square statistic is available at the \(\phi\) value being used in that M-step. Finding a new nonlinear \(\phi\) may require a separate numerical optimization.

Deriving the Random-Effect Variance Update

The part of the complete-data log likelihood involving the random effects and \(\omega^2\) is, up to terms not depending on \(\omega^2\), the sum of the Gaussian log densities for the \(N\) effects. Taking its conditional expectation replaces the sum of squared effects with its SAEM estimate \(S_{b^2}\). Thus, the part of the M-step objective that depends on \(\omega^2\) is

If \(S_{b^2}>0\), differentiation gives

The derivative is positive when \(\omega^2<S_{b^2}/N\) and negative when \(\omega^2>S_{b^2}/N\). The objective therefore increases up to that value and decreases afterward. It also tends to negative infinity as \(\omega^2\) approaches zero or infinity when \(S_{b^2}>0\). The stationary point is consequently the global maximum over positive variances.

The condition \(S_{b^2}>0\) matters. If \(S_{b^2}=0\), then \(Q_b(\omega^2)=-(N/2)\log(\omega^2)\), which tends to positive infinity as \(\omega^2\) approaches zero from above. There is no finite maximum over positive variances, and the ordinary Gaussian density is not defined at zero variance. Zero is only a limiting boundary solution in this case, not an attained maximum of the stated objective. An implementation must handle this limit explicitly, for example by imposing a strictly positive lower bound or by using a model or penalty that defines a different objective. A small positive variance floor is a computational constraint; it should not be described as the unconstrained maximum-likelihood update.

Deriving the Residual Variance Update

Under the stated independent, homoscedastic Gaussian error model, each of the \(L\) observations contributes a log-variance term and a squared residual term. Conditional expectation and SAEM approximation replace the total squared residual by \(S_{e^2}(\phi)\). Holding the prediction parameters at the value used to calculate that statistic, the variance-dependent objective is

When \(S_{e^2}(\phi)>0\), differentiation yields

This derivative is positive below \(S_{e^2}(\phi)/L\) and negative above it. Also, for a positive residual-square statistic, \(Q_e\) tends to negative infinity at both ends of the positive variance domain. Thus the stationary point is the global maximum over positive variances.

If \(S_{e^2}(\phi)=0\), then \(Q_e(\sigma^2)=-(L/2)\log(\sigma^2)\), which tends to positive infinity as \(\sigma^2\) approaches zero from above. There is no finite maximum for positive variance, and the Gaussian density is not defined at zero variance. The formula gives zero only as a limiting boundary update in this case; it is not an attained maximum under the ordinary positive-variance Gaussian model. An implementation must explicitly constrain or otherwise handle that limit. In a nondegenerate continuous-error model, an exactly zero expected residual-square statistic is unusual, but numerical degeneracy or a defective statistic can produce it.

Using Both Updates in an SAEM M-Step

When both Gaussian variance parameters are in the model, the random-effect and residual contributions can be maximized separately if the other parameters are held fixed and the relevant sufficient statistics have been accumulated. This separability follows because the random-effect variance enters the random-effect density, while the residual variance enters the conditional observation density. It does not mean that every model parameter can be updated separately: a nonlinear prediction parameter can affect the residual-square statistic.

Diagnostics, Alternatives, and Limitations

A closed-form update is simple only if the statistic and its assumptions are right. A negative accumulated squared-statistic is incompatible with its definition and signals an implementation or bookkeeping error. A variance with units inconsistent with the model, a denominator that counts the wrong units, or a residual statistic calculated using a different prediction parameter value can likewise produce a numerically plausible but invalid update. Check the definitions of the stored statistics, not just the resulting parameter traces.

Very small positive statistics yield very small positive variance updates. They may reflect data that support little variability, but they can also arise from poor latent-effect exploration, a degenerate approximation, or a coding error. The closed-form formula alone cannot distinguish these explanations. Compare the statistic over iterations and verify that it is being updated from the intended SAEM draws and gain schedule. Conversely, a noisy statistic can create a noisy parameter update; a closed-form M-step does not remove Monte Carlo variability in the E-step.

These formulas rely on specific distributional structure. If random effects have a nonzero estimated mean, if their covariance is not the scalar variance shown here, or if residual variance depends on the observation or prediction, the likelihood terms and sufficient statistics change. The scalar formulas should not be applied to those cases without deriving the corresponding objective. Likewise, if zero is reached as a limiting update, choosing a positive floor changes the feasible parameter space and should be reported as such.

The main alternative for parameters without a closed-form maximizer is a numerical M-step: optimize the expected complete-data objective over those parameters while holding the needed SAEM statistics fixed. A direct numerical optimizer may also be used for a variance block, but it should reproduce the same interior maximizer under the assumptions above. Numerical optimization does not fix a misspecified error model or an undefined boundary objective. In SAEM, the closed-form updates are convenient components of the M-step, not a guarantee that every parameter has a closed-form update or that every stochastic iteration increases the observed-data likelihood.