Gaussian Random-Effect Statistics in SAEM

Why Gaussian Random-Effect Statistics Matter

In SAEM, the latent random effects are simulated as part of the complete data, while the population parameters are updated using statistics of those simulated effects. For Gaussian random effects, the relevant information is compact: a first-moment statistic and a second-moment statistic. Their stochastic averages provide the quantities needed for the mean and covariance updates, even when the random effects themselves change from iteration to iteration.

The key is to average the sufficient statistics, not to confuse them with estimates computed from a single simulated set of random effects. The resulting covariance update also needs a qualification: its centered matrix must be positive definite for the unconstrained Gaussian likelihood to have a finite maximizer in the positive-definite covariance parameter space.

Definition: Suppose $N$ individuals have independent random-effect vectors $b_i\in\mathbb{R}^q$, with $b_i\sim N(\mu,\Omega)$. Here $\mu$ is the $q$-dimensional population mean and $\Omega$ is a symmetric positive-definite $q\times q$ covariance matrix. The Gaussian random-effect sufficient statistics are $$ T_1(b)=\sum_{i=1}^{N}b_i,\qquad T_2(b)=\sum_{i=1}^{N}b_i b_i^\mathsf{T}. $$

These statistics are sufficient for estimating $\mu$ and $\Omega$ from the complete set of random effects under this Gaussian model. In an NLME model, the complete-data likelihood also includes contributions from the observations conditional on the effects. Those contributions may involve other population parameters, but the Gaussian random-effect part depends on $\mu$ and $\Omega$ through $T_1$ and $T_2$.

From Gaussian Effects to M-Step Updates

Ignoring terms that do not depend on $\mu$ or $\Omega$, the complete-data log-likelihood for the random effects is

\ell(\mu,\Omega) =-\frac{N}{2}\log|\Omega| -\frac{1}{2}\sum_{i=1}^{N}(b_i-\mu)^\mathsf{T}\Omega^{-1}(b_i-\mu).

The mean and covariance updates can be expressed using the raw statistics. Define the centered sum-of-products matrix

C=T_2-\frac{1}{N}T_1T_1^\mathsf{T}.

Expanding the centered quadratic sum gives

\sum_{i=1}^{N}(b_i-\mu)(b_i-\mu)^\mathsf{T} =T_2-T_1\mu^\mathsf{T}-\mu T_1^\mathsf{T}+N\mu\mu^\mathsf{T} =C+N(\mu-\bar b)(\mu-\bar b)^\mathsf{T}, \qquad \bar b=\frac{T_1}{N}.

For fixed positive-definite $\Omega$, the last term contributes a nonnegative quadratic form to the negative log-likelihood, minimized at $\mu=\bar b$. At that mean, the centered sum is $C$, and the covariance-dependent part of the log-likelihood becomes

\ell(\bar b,\Omega) =-\frac{N}{2}\log|\Omega|-\frac{1}{2}\operatorname{tr}(\Omega^{-1}C)+\text{constant}.

If $C$ is positive definite, differentiating with respect to $\Omega$ gives the stationary point $\Omega=C/N$. This is the unique maximum in the positive-definite covariance space. The divisor is $N$, because this is a maximum-likelihood update, not the unbiased sample covariance formula.

Important qualification: The unconstrained Gaussian maximum-likelihood covariance update is $\widehat{\Omega}=C/N$ when $C$ is positive definite. If $C$ is singular, there is no finite maximizer over positive-definite covariance matrices: the likelihood can increase without bound as covariance eigenvalues approach zero in directions with no observed centered variation. A constrained or regularized update, or other model-specific handling, is then required.

The qualification matters in practice. Since $C$ is formed from $N$ centered vectors, its rank is at most $N-1$. Thus, if $q\geq N$, it cannot be positive definite. Singular or nearly singular statistics can also occur when the simulated random effects have little variation in one or more directions.

Worked Example: Scalar Random Effects and the Likelihood Divisor

Suppose a hypothetical complete-data set contains $N=3$ scalar random effects, measured in arbitrary effect units: $b_1=0.2$, $b_2=-0.1$, and $b_3=0.5$. Calculate the Gaussian maximum-likelihood mean and variance updates.

The first-moment statistic is $T_1=0.2-0.1+0.5=0.6$. The second-moment statistic is $T_2=0.2^2+(-0.1)^2+0.5^2=0.04+0.01+0.25=0.30$. Thus the estimated mean is $ \widehat{\mu}=T_1/N=0.6/3=0.20 $ effect units.

The centered sum of squares is $C=T_2-T_1^2/N=0.30-0.6^2/3=0.30-0.12=0.18$ squared effect units. The maximum-likelihood variance update is $ \widehat{\Omega}=C/N=0.18/3=0.0600 $ squared effect units.

For comparison, the unbiased sample variance would be $0.18/(3-1)=0.0900$ squared effect units. That is a different estimator, with a different objective. The Gaussian likelihood update uses divisor $N$, as derived above.

Stochastic Approximation of the Statistics

At SAEM iteration $k$, let $b_i^{(k)}$ denote the simulated random effect for individual $i$. The statistics calculated from that iteration’s simulated effects are

T_1^{(k)}=\sum_{i=1}^{N}b_i^{(k)},\qquad T_2^{(k)}=\sum_{i=1}^{N}b_i^{(k)}b_i^{(k)\mathsf{T}}.

SAEM carries running statistics forward by a stochastic approximation step. With gain $\gamma_k$, the updates are

S_{1,k}=(1-\gamma_k)S_{1,k-1}+\gamma_k T_1^{(k)},\qquad S_{2,k}=(1-\gamma_k)S_{2,k-1}+\gamma_k T_2^{(k)}.

These are weighted running statistics, not necessarily the raw statistics of the current draw. The corresponding population updates are

\widehat{\mu}_k=\frac{S_{1,k}}{N},\qquad C_k=S_{2,k}-\frac{1}{N}S_{1,k}S_{1,k}^\mathsf{T},\qquad \widehat{\Omega}_k=\frac{C_k}{N},

provided $C_k$ is positive definite for the unconstrained covariance maximum. Equivalently, define normalized running moments $m_{1,k}=S_{1,k}/N$ and $m_{2,k}=S_{2,k}/N$; then $\widehat{\mu}_k=m_{1,k}$ and the covariance candidate is $m_{2,k}-m_{1,k}m_{1,k}^\mathsf{T}$. Both formulations give the same update when normalization is applied consistently.

The distinction between the second moment and the covariance is essential. $S_{2,k}/N$ averages outer products about zero; it is not itself the covariance unless the mean is zero and fixed. The covariance requires subtracting the outer product of the estimated mean. This centering is performed after forming the running moments.

Worked Example: A Two-Dimensional Covariance Update

Consider three hypothetical simulated random-effect vectors, each with two components measured in effect units: $b_1=(1,0)^\mathsf{T}$, $b_2=(0,1)^\mathsf{T}$, and $b_3=(1,1)^\mathsf{T}$. First calculate their mean and covariance candidate.

Their sum is $T_1=(2,2)^\mathsf{T}$, so $\bar b=(2/3,2/3)^\mathsf{T}$. Their outer products sum to $T_2=\begin{pmatrix}2&1\\1&2\end{pmatrix}$. Therefore

\frac{T_2}{3} -\bar b\bar b^\mathsf{T} = \begin{pmatrix}2/3&1/3\\1/3&2/3\end{pmatrix} - \begin{pmatrix}4/9&4/9\\4/9&4/9\end{pmatrix} = \begin{pmatrix}2/9&-1/9\\-1/9&2/9\end{pmatrix}.

This matrix is positive definite: its eigenvalues are $1/9$ and $1/3$, both positive. Thus the Gaussian maximum-likelihood covariance update is the displayed matrix. Its diagonal entries are the estimated component variances, $2/9$ squared effect units each, and its off-diagonal entries are the estimated covariance, $-1/9$ squared effect units.

This check is not optional when applying the unconstrained update. A symmetric matrix is not automatically a valid covariance matrix; it must also be positive definite in the parameter space specified here.

Worked Example: One Stochastic-Averaging Step

For a hypothetical scalar random effect with $N=2$ individuals, suppose the normalized previous running moments are $m_{1,k-1}=0.10$ effect units and $m_{2,k-1}=0.05$ squared effect units. At the current iteration the simulated effects are $-0.2$ and $0.6$, and the gain is $\gamma_k=0.25$. Find the updated moments and covariance candidate.

The current draw gives a first moment of $(-0.2+0.6)/2=0.20$ effect units and a second moment of $((-0.2)^2+0.6^2)/2=(0.04+0.36)/2=0.20$ squared effect units. Applying the same gain to both statistics gives

m_{1,k}=0.75(0.10)+0.25(0.20)=0.125,\qquad m_{2,k}=0.75(0.05)+0.25(0.20)=0.0875.

The mean update is $0.125$ effect units. The covariance candidate is $m_{2,k}-m_{1,k}^2=0.0875-0.125^2=0.0875-0.015625=0.071875$ squared effect units. It is positive, so this scalar covariance update is a valid finite maximum-likelihood update for these running statistics.

Checks, Failure Modes, and Alternatives

Several checks help identify an incorrect or unusable update:

Check dimensions. $S_{1,k}$ must be a $q$-vector and $S_{2,k}$ a $q\times q$ matrix. The covariance has squared effect units and the same units as an outer product.
Check symmetry. $S_{2,k}$ and $C_k$ should be symmetric. Small numerical asymmetry can arise from computation; substantial asymmetry indicates an implementation or data-handling problem.
Check positive definiteness. A zero or negative eigenvalue means the unconstrained positive-definite covariance update is not available from that centered statistic. A tiny positive eigenvalue signals a poorly conditioned estimate, even though it is technically positive definite.
Use matching moments. Center the second moment with the mean derived from the same running first-moment statistic. Mixing statistics from different iterations or normalizations changes the covariance calculation.

The raw-moment representation $T_1,T_2$ is compact and convenient for stochastic averaging. An alternative is to average centered quantities around a separately updated mean, or to work with a parameterization that enforces covariance constraints. These alternatives can be useful computationally, but they must preserve the intended Gaussian likelihood update; they are not interchangeable without checking how their averaging and centering are defined. In particular, replacing the divisor $N$ by $N-1$ targets an unbiased sample covariance, not the Gaussian maximum-likelihood update.

Common pitfall: “The centered sum-of-products divided by $N$ is always the covariance MLE because it solves the stationary equation.” The stationary equation gives a valid finite maximum only when the centered matrix is positive definite. If it is singular, the likelihood has no finite maximizer over positive-definite covariance matrices; use an appropriate constraint, regularization, or model-specific treatment rather than reporting a singular matrix as an unconstrained MLE.

These statistics matter because they separate the latent-effect simulation step from the population-parameter calculation: SAEM can average low-dimensional moments rather than retain the full history of every simulated effect. Their limitation is equally important: stochastic averaging does not guarantee that every finite-iteration covariance candidate is well-conditioned or even positive definite. The candidate must still be checked before it is interpreted as a valid unconstrained Gaussian maximum.

Check Your Understanding

Answer each question, then compare with the solutions below.

For $N$ scalar Gaussian random effects, write the first- and second-moment statistics and the covariance update in terms of them.
Why does the Gaussian maximum-likelihood covariance use divisor $N$, rather than $N-1$?
For two-dimensional effects and $N=2$, why must the centered sum-of-products matrix be singular when the mean is estimated?
A running covariance candidate is symmetric but has one zero eigenvalue. Can it be reported as a finite unconstrained maximum-likelihood covariance in the positive-definite parameter space?

Solutions. (1) The statistics are $T_1=\sum_i b_i$ and $T_2=\sum_i b_i^2$; the mean is $T_1/N$, and the variance candidate is $(T_2-T_1^2/N)/N$, if positive. (2) Maximizing the Gaussian likelihood gives the centered sum of squares divided by the number of observations, $N$; $N-1$ is the divisor for the unbiased sample-variance estimator instead. (3) Centering $N=2$ vectors gives two deviations whose sum is zero, so they span at most one direction; the centered matrix has rank at most $N-1=1$, less than dimension $q=2$. (4) No. The matrix is singular, so there is no finite maximizer over positive-definite covariance matrices; an appropriate constrained or regularized treatment is needed.

Gaussian Random-Effect Statistics in SAEM

What You'll Learn