Residual-Error Updates From SAEM Statistics

Why the Residual-Variance Update Needs Its Own Statistic

In a nonlinear mixed-effects model, the prediction for an observation depends on the subject’s latent random effects. SAEM draws those effects while approximating expectations in the complete-data likelihood. A residual-variance update therefore needs the squared discrepancy between each observation and its model prediction, evaluated conditional on the latent-effect draw. It is not an update of the random-effect covariance: that describes variation between subjects in their latent effects, whereas residual variance describes observation-level variation around the conditional prediction.

The distinction matters because both sources of variation can make observed values differ across subjects. If their contributions are mixed up, the residual-error parameter can be too large or too small, and the model’s uncertainty for individual observations will be misrepresented.

Conditional Error Models and Their Assumptions

For an additive-error model, conditional on the random effects, observations have constant residual variance in the response’s units:

This model assumes conditionally independent observation errors with mean zero and variance \(\sigma_a^2\), given the random effects and model parameters. The additive residual standard deviation \(\sigma_a\) has the same units as \(y_{ij}\).

For a proportional-error model, write the observation as \(y_{ij}=f_{ij}(\theta,b_i)(1+\epsilon_{ij})\), with \(\epsilon_{ij}\sim N(0,\sigma_p^2)\). Equivalently, conditional on the random effects:

Here the residual standard deviation is \(\sigma_p|f_{ij}|\), so \(\sigma_p\) is dimensionless and is often interpreted as a relative residual scale. This formulation assumes conditionally independent, mean-zero proportional errors. The prediction used as the scale must be nonzero for the relative-residual expression below to be defined. When predictions can be zero or very close to zero, a purely proportional model may be unsuitable.

Both error models condition on \(b_i\). The random effects have their own population distribution and covariance parameters. In the residual update, a draw of \(b_i\) is used to calculate \(f_{ij}\); the draw itself is not a residual and is not added to the residual sum of squares.

Deriving the Closed-Form Updates

Consider the expected complete-data log likelihood for the residual parameter, holding the structural parameters and the random-effect distribution parameters fixed. “Complete data” here means the observed \(y_{ij}\) together with the latent \(b_i\). For an additive model, the terms involving the residual variance \(v_a=\sigma_a^2\) are

where \(C\) does not depend on \(v_a\). Write \(A_a\) for the conditional expectation of the summed squared residuals. Differentiating with respect to \(v_a\), setting the derivative to zero, and solving gives the update:

The same reasoning applies to proportional error, but each squared residual is divided by the squared prediction. For fixed predictions, set \(A_p=\operatorname{E}[\sum_{i,j}e_{ij}^2/f_{ij}^2\mid y]\). The variance-dependent terms in the log likelihood are \(-N\log(v_p)/2-A_p/(2v_p)\), where \(v_p=\sigma_p^2\). Thus:

The denominators are the same \(N\) because each observation contributes one Gaussian residual term. These are variance updates, not standard-deviation updates: compute the mean squared residual on the relevant scale first, then take its square root if a standard deviation is required. In an exact EM M-step, the expectation is under the conditional distribution of the random effects given the data at the preceding parameter values. SAEM approximates this expectation through latent-effect simulation and stochastic approximation.

Putting the Residual Statistic Into SAEM

At each SAEM iteration, latent random effects are simulated from an appropriate conditional distribution, and the resulting complete-data statistic is incorporated using a stochastic-approximation gain. For a single simulated draw of all subject effects at iteration \(k\), define the additive and proportional residual sums:

The residual statistic is evaluated using the model predictions at the structural parameter values used for that residual-statistic calculation. If structural parameters are being updated as well, the implementation must evaluate the residual statistic consistently with its M-step; it must not treat a statistic computed for one set of predictions as if it automatically applied to another. With fixed \(N\), the Robbins–Monro updates for the running sums are:

Here \(\gamma_k\) is the gain at iteration \(k\), and the statistics are initialized and updated according to the SAEM schedule. The variance updates use \(A_{a,k}/N\) or \(A_{p,k}/N\), respectively. The fixed count \(N\) is not itself being estimated or averaged. If multiple conditional draws are used to estimate the statistic at an iteration, their contribution should be averaged consistently before applying the gain.

What the Update Does—and Does Not—Measure

A residual statistic is built from observation-level discrepancies after the conditional prediction has accounted for the simulated subject effects. For example, suppose two subjects have different predicted concentrations because their simulated random effects differ. The difference between those predictions is part of the model’s between-subject structure. Only the remaining observed-minus-predicted differences contribute to the conditional residual sum.

This separation also explains why the number of observations, rather than the number of subjects, appears in the denominator. Every conditionally independent observation contributes a residual likelihood term. Repeated observations from one subject share that subject’s latent effects, but they still each contribute their own conditional residual term under the stated model.

A useful implementation check is to inspect the quantities entering the update. Additive residuals should be in response units, and their squares in squared response units. Proportional residuals should be dimensionless ratios, and their squares dimensionless. A units mismatch often signals that predictions, residuals, or weights have been combined incorrectly.

Diagnostics, Alternatives, and Limitations

The closed-form update is valid under its specified conditional Gaussian error model and for fixed predictions during the residual-variance maximization. It does not establish that the chosen error model fits the data. Examine conditional residuals against predictions, time, and other relevant covariates. Under an additive model, the residual spread should be roughly stable across the prediction range. Under a proportional model, the relative residual spread should be roughly stable; absolute residual spread may increase with prediction magnitude.

Patterns in these plots can indicate a misspecified error structure, a systematic trend in the mean model, outliers, or dependence not represented by the conditional independence assumption. A very small proportional prediction deserves particular attention: dividing by it can give a single observation extreme influence. A zero prediction makes the proportional statistic undefined. Do not silently drop such observations or add an arbitrary small constant without defining and justifying a different error model.

Additive error is a natural choice when the absolute residual scale is stable, while proportional error is useful when relative variability is more nearly stable. A combined additive-plus-proportional model can represent both a nonzero baseline error scale and variability that grows with the prediction, but it has a different variance formula and requires its own update. It is not obtained by applying both formulas above independently to the same residuals. Alternatively, if the conditional variance has no closed-form update under a chosen model, its M-step may require numerical optimization of the expected complete-data log likelihood.

SAEM adds Monte Carlo variability to the statistic at finite iterations. A noisy variance trace may reflect the conditional simulation rather than a change in the underlying fit, so interpret the update together with the gain schedule and the residual-statistic trace. The update also relies on the model prediction being computed correctly for each latent-effect draw. It cannot compensate for an incorrectly specified mean function or for a mismatch between the implemented likelihood and the assumed error model.

• Common mistake: “The random-effect draws vary across subjects, so I should add their sample variance to the residual sum of squares.”
• Correction: Random-effect variation belongs to the distribution of \(b_i\). For each draw, calculate the conditional prediction \(f_{ij}(\theta,b_i)\), then form \(y_{ij}-f_{ij}(\theta,b_i)\). Only those conditional residuals enter the residual-error update.