Choosing Sufficient Statistics for SAEM

Why the Choice of Statistics Matters

In SAEM, the latent random effects are simulated or otherwise approximated as part of an iterative calculation. The algorithm does not need to retain every detail of every latent-effect draw if the complete-data likelihood depends on those draws through a smaller collection of summaries. Choosing the summaries correctly is therefore not just a matter of computational convenience: omitting a term that the likelihood needs can change the M-step and the resulting estimates.

The key distinction is between a statistic that summarizes the complete data for the likelihood and a quantity that is merely useful to monitor an iteration. For an NLME model, the complete data consist of the observed responses together with the latent random effects. The observed responses are fixed during fitting; SAEM approximates conditional expectations of the parts of the complete-data likelihood that involve latent effects.

Here \(y_i\) denotes the response vector for subject \(i\), and \(b_i\) denotes that subject’s latent random-effect vector. Let \(N\) be the number of subjects and \(n_i\) the number of observations for subject \(i\). The examples use Gaussian distributions with independent subjects and, where stated, conditionally independent residuals. These assumptions matter: a covariance parameter for correlated errors, for example, generally needs different summaries from a single residual variance.

Separate the Random-Effect and Residual Contributions

Suppose the random effects have a multivariate Gaussian distribution with mean zero and covariance matrix \(\Omega\), and the responses are conditionally Gaussian given the random effects:

Assume subjects are independent, observations are conditionally independent given their subject’s random effects, and \(\sigma^2\) is a common residual variance. The complete-data log-likelihood, omitting constants that do not depend on parameters, is

This expression shows two distinct sources of summaries. The Gaussian random-effect density depends on the latent effects through their quadratic forms; the residual density depends on them through the residual sum of squares. If \(\Omega\) is unrestricted and \(\sigma^2\) is common, useful quantities include \(\sum_i b_i b_i^\mathsf{T}\) and the residual sum of squares, evaluated at the relevant population parameters. If the mean of the random effects is also estimated, \(\sum_i b_i\) is needed as well as the quadratic summary.

The Gaussian assumption alone does not determine one universal list of SAEM statistics. The covariance structure, residual variance structure, and form of \(f_{ij}\) all affect what the likelihood requires. In particular, the residual sum of squares is a compact statistic for estimating a common variance when the mean parameters are fixed, but it is not automatically a fixed, parameter-free summary that suffices for every nonlinear mean parameter.

A Short Derivation: Why Gaussian Variances Need Quadratic Summaries

For independent zero-mean Gaussian random effects with common covariance \(\Omega\), the density of \(b_i\) is proportional to

Multiplying the subject densities adds their log-likelihood contributions. The quadratic terms can be rewritten using the trace identity \(b_i^\mathsf{T}\Omega^{-1}b_i=\operatorname{tr}(\Omega^{-1}b_i b_i^\mathsf{T})\):

Thus, for an unrestricted covariance matrix, the complete-data likelihood uses the latent effects through the matrix sum \(\sum_i b_i b_i^\mathsf{T}\). If the covariance is constrained, fewer distinct components may be needed; if the covariance is diagonal, for instance, the likelihood uses the sums of squared effects component by component. This follows from the likelihood itself, rather than from a general rule that every model must store the entire matrix.

Residual Structure Determines the Residual Summary

For independent residuals with known positive variance multipliers \(v_{ij}\), suppose \(\varepsilon_{ij}\mid b_i\sim N(0,\sigma^2v_{ij})\). The residual log-likelihood contains both a weighted squared-residual term and a determinant term:

If the \(v_{ij}\) values are known and fixed, the parameter-dependent residual summary for \(\sigma^2\) is the weighted residual sum of squares. The unweighted sum is not interchangeable with it. The log-\(v_{ij}\) terms are constant with respect to \(\sigma^2\), but must be retained if the variance multipliers themselves are parameters.

The same principle applies to other Gaussian residual structures. When residuals are correlated, the quadratic term is a covariance-weighted form, such as \(r_i^\mathsf{T}R_i^{-1}r_i\), and the likelihood also includes \(\log|R_i|\). A scalar unweighted residual sum of squares generally cannot replace those terms. The covariance model determines which summaries and parameter-dependent calculations are needed.

Nonlinear Means: A Summary May Depend on the Parameters

When \(f_{ij}(b_i;\theta)\) is nonlinear, the residual sum of squares is a function of both the latent effects and the population parameters. For a fixed \(\theta\), the sum of squared residuals is the quantity entering the likelihood for a common Gaussian residual variance. But if SAEM must estimate \(\theta\) too, there may be no finite, parameter-free collection of simple statistics that reduces the complete-data likelihood to a standard exponential-family form.

For example, the expansion

shows why. The response-only term \(\sum y_{ij}^2\) is fixed, but the other terms depend on the nonlinear function and its parameters. The conditional expectation needed for an M-step may therefore require evaluating the complete-data log-likelihood over candidate parameter values, or approximating a suitable function of the latent effects, rather than accumulating only one universal residual statistic.

Practical Selection, Diagnostics, and Limitations

A reliable way to choose SAEM summaries is to start from the complete-data log-likelihood and identify the terms involving each parameter block. Then determine exactly which functions of the latent effects occur in those terms. This avoids assuming that a familiar Gaussian model always calls for the same statistic.

Typical warning signs are a variance update that uses an unweighted residual sum where the model specifies weights, a covariance update that keeps only diagonal terms when the covariance is unrestricted, or an M-step that cannot reconstruct the complete-data objective from the accumulated quantities. A further diagnostic is to compare the full complete-data log-likelihood with its proposed summary-based representation at several parameter values. Agreement at one point alone is not enough if parameter-dependent terms were accidentally omitted.

An alternative to storing a compact set of sufficient statistics is to retain or recompute richer information, such as complete-data likelihood contributions at candidate values during a numerical M-step. This can be necessary for nonlinear mean structures, though it may cost more computation or storage. The compact-statistic approach is most direct when the complete-data likelihood has a tractable exponential-family form. It is not a license to compress away terms simply because a model is Gaussian.

These rules assume the stated Gaussian distributions and the modeled independence or covariance structure. A model with unknown residual correlation, unknown variance multipliers, or a non-Gaussian error distribution needs summaries derived from its own likelihood. Likewise, a random-effect covariance constraint changes which components of the quadratic summary are needed. The correct summaries follow from the specified model, not from the name of the algorithm.