Why update the covariance from cross-products?
In a nonlinear mixed-effects model, each individual has a vector of random effects describing deviations from the population parameters. For example, one component might describe an individual’s deviation in log clearance and another their deviation in log volume. The population covariance matrix describes how these deviations vary across individuals and whether they tend to vary together. In SAEM, the random effects are latent, so this matrix is not updated by calculating an ordinary covariance from observed individual effects. Instead, its update uses conditional random-effect cross-products, estimated from the simulated draws used by the algorithm.
The focus here is the covariance M-step: why each off-diagonal entry is an average cross-product, how SAEM accumulates that information, and what the resulting entry does and does not mean. The introductory material on Gaussian random-effect statistics in SAEM establishes the relevant sufficient-statistic framework; the derivation below applies it specifically to a full covariance matrix.
The zero mean is part of this model: systematic population parameters are represented separately, while \(\eta_i\) describes deviations from them. Consequently, the covariance is based on second moments \(E(\eta_i\eta_i^T)\), not on subtracting an estimated mean of the latent effects. If a model instead estimates a nonzero random-effect mean, its covariance update must account for that mean; that is not the model considered here.
Deriving the covariance M-step
The complete-data contribution of the random effects to the log likelihood, omitting constants that do not involve \(\Omega\), is
In the EM M-step, replace the latent quantities by their conditional expectations given the observations and the current parameter values. Define the average conditional second-moment matrix
The part of the expected complete-data log likelihood that depends on \(\Omega\) becomes
Here the trace identity \(\eta_i^T\Omega^{-1}\eta_i=\operatorname{tr}(\Omega^{-1}\eta_i\eta_i^T)\) turns the quadratic form into a matrix cross-product. Differentiating with respect to the symmetric matrix \(\Omega\), the stationary condition is
Multiplying on the left and right by \(\Omega\) gives \(\Omega=A\). When \(A\) is positive definite, this is the maximizer over positive-definite covariance matrices. Thus the exact EM update is
This result also explains the SAEM update. With current simulated draws \(\eta_i^{(k)}\), a one-draw estimate of the complete-data statistic is \(T_k=N^{-1}\sum_i\eta_i^{(k)}\eta_i^{(k)T}\). If \(L\) draws per individual are used at iteration \(k\), average the outer products over those draws to form \(T_k\). The stochastic-approximation recursion is \(S_k=S_{k-1}+\gamma_k(T_k-S_{k-1})\), and the covariance update is \(\Omega_k=S_k\). The gain \(\gamma_k\) controls how much the current simulated statistic changes the running statistic.
Worked covariance updates
Worked Example: Updating a Two-Effect Covariance From Simulated Draws
Suppose two hypothetical individuals have two simulated draws each for a pair of random effects, \(\eta_i=(\eta_{i1},\eta_{i2})^T\). All values are illustrative and are expressed as deviations on a log scale. For individual 1 the draws are \((1,1)^T\) and \((0,0)^T\); for individual 2 they are \((2,1)^T\) and \((0,1)^T\). Estimate each individual’s conditional second moment by averaging the two outer products, then average across the two individuals.
For individual 1, the outer products are \(\begin{pmatrix}1&1\\1&1\end{pmatrix}\) and \(\begin{pmatrix}0&0\\0&0\end{pmatrix}\). Their average is \(\begin{pmatrix}0.5&0.5\\0.5&0.5\end{pmatrix}\). For individual 2, the outer products are \(\begin{pmatrix}4&2\\2&1\end{pmatrix}\) and \(\begin{pmatrix}0&0\\0&1\end{pmatrix}\); their average is \(\begin{pmatrix}2&1\\1&1\end{pmatrix}\).
The covariance update is therefore
The diagonal entries are the estimated variances of the two log-scale random effects. The off-diagonal covariance is \(0.75\), in squared log-deviation units. The corresponding correlation is \(0.75/\sqrt{1.25(0.75)}\approx0.7746\), which is dimensionless. The matrix is positive definite because its determinant is \(1.25(0.75)-0.75^2=0.375>0\), and its first diagonal entry is positive.
Worked Example: Why the Product of Conditional Means Is Not Enough
For one hypothetical individual, suppose the simulated conditional draws for two random effects are \((1,1)^T\) and \((-1,-1)^T\), each with equal weight. The average first component is \(0\), as is the average second component, so multiplying the two conditional means gives \(0\times0=0\).
The cross-product statistic instead averages the products within each draw:
The same distinction appears in the full second-moment matrix:
Thus, zero conditional means do not imply a zero conditional cross-product. The two effects move together in these draws, and that joint information contributes to the off-diagonal covariance update. Using the outer product of the conditional-mean vector would discard this variation around the means.
Worked Example: A SAEM Covariance-Statistic Step
Suppose a hypothetical SAEM iteration starts with running covariance statistic \(S_{k-1}=\begin{pmatrix}1&0.2\\0.2&0.5\end{pmatrix}\). The current simulated complete-data cross-product statistic, already averaged over individuals, is \(T_k=\begin{pmatrix}4&1.2\\1.2&1\end{pmatrix}\), and the gain is \(\gamma_k=0.25\). Apply the stochastic-approximation recursion.
The old statistic receives weight \(1-\gamma_k=0.75\) and the current statistic receives weight \(0.25\):
The updated covariance estimate is \(\Omega_k=S_k\). Its off-diagonal covariance is \(0.45\), and its correlation is \(0.45/\sqrt{1.75(0.625)}\approx0.4303\). The determinant is \(1.75(0.625)-0.45^2=0.89125>0\), so this particular two-dimensional update is positive definite. The calculation is a weighted update of the sufficient statistic, not a claim that the current iteration alone estimates the covariance as \(S_k\).
Interpreting dependence and checking the update
For random effects \(\eta_{i1}\) and \(\eta_{i2}\), the population covariance is \(\Omega_{12}=E(\eta_{i1}\eta_{i2})\) under the zero-mean model. A positive value indicates that individuals with above-population deviations in one component tend, in the population model, to have above-population deviations in the other. A negative value indicates a tendency for opposite-signed deviations. The magnitude depends on the scales of the two components, so correlation is often more useful for describing strength:
This interpretation concerns the joint distribution of individual random effects. It is not by itself a claim that observed outcomes are correlated in the same way, nor does it establish that changing one individual parameter causes another to change. The nonlinear structural model and residual variation also affect observed outcomes.
Several checks help identify a faulty or unstable update:
- Symmetry: an outer product is symmetric, so the accumulated statistic should also be symmetric. A material mismatch between \(\Omega_{ab}\) and \(\Omega_{ba}\) suggests an implementation or numerical error.
- Diagonal entries: estimated variances must be nonnegative. A negative diagonal entry is incompatible with a covariance matrix and signals an error.
- Positive definiteness: a covariance used as a nonsingular Gaussian covariance must be positive definite. Finite simulated draws can yield a singular cross-product matrix if they do not span the full random-effect space. Numerical roundoff can also cause problems. The earlier material on maintaining positive-definite random-effect covariance addresses this constraint; symmetry alone is not enough.
- Scale and stability: an off-diagonal estimate can be hard to interpret when random-effect scales differ greatly. Inspect the implied correlation as well as the covariance, and check whether the estimate stabilizes rather than fluctuating sharply with simulated draws.
- Simulation quality: the update is only as representative as the conditional draws. Poor exploration can omit parts of the conditional distribution and distort its cross-products. A noisy off-diagonal trace is not automatically evidence of strong population dependence.
Alternatives, assumptions, and limitations
The closed-form update applies when the random effects are Gaussian with zero mean and a shared, unrestricted covariance matrix, and when the conditional distribution used by EM or SAEM is represented adequately. The ordinary EM formula uses exact conditional expectations; SAEM substitutes a stochastic approximation to the sufficient statistic. If the covariance is constrained—for example, required to be diagonal or to follow a lower-dimensional structure—the unrestricted formula is not automatically the appropriate constrained maximizer.
| Approach | What it uses | Main trade-off |
|---|---|---|
| Full covariance M-step | All conditional second moments and cross-products | Represents dependence, but estimates more entries and may be unstable with limited information |
| Diagonal covariance model | Second moments for each component; off-diagonal entries are fixed at zero | Simpler, but assumes no population covariance between distinct random effects |
| Covariance of conditional means | One posterior mean vector per individual | Easy to calculate, but omits conditional uncertainty and is generally not the EM covariance update |
The full update is preferable only when the model and data can support the additional covariance parameters. A small or erratic off-diagonal estimate may reflect weak information, limited conditional simulation, or a genuinely weak association; the update alone cannot distinguish those explanations. Likewise, a stable covariance estimate is not proof that the Gaussian random-effect distribution or the assumed model structure is adequate. Interpretation should remain tied to the model and its diagnostics.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- For a zero-mean Gaussian random-effect model, what quantity contributes to the covariance M-step for individual \(i\): the conditional mean vector or the conditional expected outer product?
- A two-effect covariance estimate has variances \(4\) and \(9\), and covariance \(-3\). Calculate the correlation and interpret its sign.
- A simulated-draw average gives conditional means \((0,0)^T\), but the average cross-product for the two components is \(0.6\). Which value enters the covariance update?
- With gain \(\gamma_k=0.2\), the previous covariance statistic is \(\begin{pmatrix}2&0.4\\0.4&1\end{pmatrix}\) and the current statistic is \(\begin{pmatrix}3&-0.1\\-0.1&2\end{pmatrix}\). Calculate the next statistic.
Solutions. (1) Use \(E(\eta_i\eta_i^T\mid y_i)\), the conditional expected outer product. (2) The correlation is \(-3/\sqrt{4(9)}=-3/6=-0.5\); the random effects have a negative population association in the model, with correlation measured on a unit-free scale. (3) The value \(0.6\) enters as the off-diagonal second-moment contribution; multiplying the conditional means would incorrectly give zero. (4) Apply \(S_k=0.8S_{k-1}+0.2T_k\), giving \(S_k=\begin{pmatrix}1.6&0.32\\0.32&0.8\end{pmatrix}+\begin{pmatrix}0.6&-0.02\\-0.02&0.4\end{pmatrix}=\begin{pmatrix}2.2&0.30\\0.30&1.2\end{pmatrix}\).