A Worked EM Update for a Gaussian Random-Intercept Model

Why Work Through a Complete EM Iteration?

In a nonlinear mixed-effects model, the individual random effects are unobserved. The EM algorithm handles this by taking conditional expectations of complete-data quantities, then maximizing the resulting expected complete-data log-likelihood. A Gaussian random-intercept model is a useful setting in which to see each calculation explicitly: the conditional distribution and the M-step updates both have closed forms.

The calculations below use maximum likelihood, not an unbiased variance estimator. The E-step expectations are taken under the current parameter values; the M-step then updates the population parameters using those expectations. Keeping those two roles distinct is essential to getting a correct iteration.

Model and Assumptions

The fixed intercept is \(\beta\), measured in the units of \(y\); the random intercept \(b_i\) is a subject-specific deviation in the same units. The zero mean of \(b_i\) defines the population intercept. The assumptions include a common residual variance and a common random-intercept variance across subjects. Conditional on the random effects, the observations are independent Gaussian measurements.

At iteration \(t\), write the current parameter values as \(\theta^{(t)}=(\beta^{(t)},\tau^{2(t)},\sigma^{2(t)})\). The E-step conditions on the observed responses and these current values. The formulas below rely on the stated Gaussian model; they are not the corresponding conditional moments for a general nonlinear mixed model.

Conditional Random-Effect Moments

For subject \(i\), let \(r_{ij}^{(t)}=y_{ij}-\beta^{(t)}\). The likelihood contribution involving \(b_i\), multiplied by its Gaussian prior density, is proportional to $$ \exp\left\{-\frac{1}{2}\left[\frac{b_i^2}{\tau^{2(t)}}+ \frac{\sum_{j=1}^{n_i}(r_{ij}^{(t)}-b_i)^2}{\sigma^{2(t)}}\right]\right\}. $$ Collecting the terms quadratic and linear in \(b_i\) gives conditional precision \(1/\tau^{2(t)}+n_i/\sigma^{2(t)}\). Completing the square therefore gives

Here \(\mathbf y_i\) denotes all observations from subject \(i\). The conditional second moment is not just the square of the conditional mean:

The variance term matters because the random effect remains uncertain after observing the data. In the M-step, the random-intercept variance uses the average of these conditional second moments, not the average of the squared conditional means.

Deriving the M-Step Updates

The complete-data log-likelihood, up to terms not involving the parameters, is

In the M-step, take its conditional expectation using the E-step distribution, then maximize with respect to \(\beta,\tau^2,\sigma^2\). The terms for \(\beta\) give a weighted sum of expected squared residuals. Differentiating that sum with respect to \(\beta\) sets the sum of expected residuals to zero:

For the variance terms, maximizing \(-k\log(q)/2-S/(2q)\) over \(q>0\) gives \(q=S/k\): its derivative is \(-k/(2q)+S/(2q^2)\), which vanishes at \(q=S/k\), and the expression is maximized there. Applying this result to each variance term yields the updates:

The residual update follows from \(E[(y_{ij}-\beta-b_i)^2\mid\mathbf y_i,\theta^{(t)}]=(y_{ij}-\beta-m_i^{(t)})^2+v_i^{(t)}\). The old conditional moments appear throughout these updates; only the fixed intercept inside the residual calculation is the newly updated value. The denominators are \(m\) for the random-effect variance and \(N\) for the residual variance because these are maximum-likelihood updates.

Checks, Diagnostics, and Limitations

Before interpreting an iteration, check that the conditional variances are positive, that \(E[b_i^2\mid\mathbf y_i,\theta^{(t)}]\geq (m_i^{(t)})^2\), and that the variance updates use the correct denominators. Also check units: \(\beta\) and \(m_i^{(t)}\) have response units, while \(\tau^2\), \(\sigma^2\), and the conditional second moments have squared response units.

A negative updated variance, or a second moment smaller than the squared mean, signals an arithmetic or implementation error under this model. A variance approaching zero is different: it may be a boundary estimate or a sign that the data provide little information to distinguish variation at that level. It warrants examination of the fit and assumptions, rather than silently replacing the estimate with an arbitrary positive number.

The closed forms depend on the linear Gaussian random-intercept structure. With a nonlinear observation model, the conditional distribution of random effects is generally not Gaussian, and these exact conditional moments and closed-form updates need not apply. Direct marginal-likelihood optimization is an alternative way to estimate parameters, but it requires evaluating or approximating the integral over random effects. The earlier material on the EM and SAEM algorithms explains why exact EM can be difficult in nonlinear mixed models; the present calculations are an exact, tractable special case.

A common mistake is to update \(\tau^2\) using only \(\sum_i(m_i^{(t)})^2/m\). That treats each conditional mean as if it were the known random effect. The correct update adds the uncertainty for every subject: \(\tau^{2(t+1)}=\sum_i[(m_i^{(t)})^2+v_i^{(t)}]/m\). Another mistake is to calculate the conditional moments using the newly updated parameters partway through the E-step. In an EM iteration, all conditional moments for that iteration use \(\theta^{(t)}\); the resulting expectations are then used in the M-step.

Under the stated model and an exact E-step and M-step, the observed-data likelihood does not decrease from one EM iteration to the next. That property does not guarantee rapid convergence or prove that the fitted model is appropriate. In practice, inspect the likelihood trace and parameter changes, and assess whether the assumptions of a common Gaussian random-intercept distribution and independent Gaussian residuals are plausible for the application.