What the E-Step Calculates
In an NLME model, the random effects are unobserved, but the complete-data log-likelihood treats them as if they were observed along with the response data. The E-step does not replace each subject’s random effects by a single best-fitting value. Instead, it averages the complete-data log-likelihood—or the statistics needed to represent it—over the conditional distribution of those effects given that subject’s observations and the current population-parameter estimates.
Let \(y_i\) denote subject \(i\)’s observed data, \(b_i\) that subject’s latent random-effects vector, and \(\theta\) the population parameters. Suppose subjects are independent conditional on the model parameters, with random effects drawn from their population distribution and responses generated according to the subject-specific model. Write the complete-data log-likelihood as
The current parameter value is written \(\theta^{(k)}\). The conditional distribution needed for subject \(i\) is
The denominator is the subject’s marginal density under the current parameters. The expectation in the E-step is with respect to this conditional distribution, not the unconditional population distribution of \(b_i\).
From the Complete-Data Log-Likelihood to Statistics
The earlier material on constructing the complete-data log-likelihood identifies its response and random-effect terms. The next step is to determine which parts of those terms depend on the latent effects. For example, suppose the complete-data log-likelihood can be written in exponential-family form as
where \(T(y,b)\) is a vector of complete-data statistics, \(\eta(\theta)\) is a parameter-dependent coefficient vector, and \(c(y)\) does not depend on \(\theta\). Taking the conditional expectation gives
Thus, in this case, the E-step needs the conditional expectations of the components of \(T\). The observed data \(y\) are fixed in this expectation; the uncertainty being averaged over is the uncertainty about the latent effects. The candidate parameter \(\theta\) in \(Q\) and the current parameter \(\theta^{(k)}\) in the conditional distribution have different roles.
A nonlinear mixed model does not always have a compact, parameter-independent sufficient-statistic vector of this form. A nonlinear response function can make the complete-data log-likelihood depend on \(b_i\) through functions whose form changes with \(\theta\). In that situation, the E-step still has a precise definition: calculate the conditional expectation of the complete-data log-likelihood, or of the parameter-dependent terms needed to evaluate it. The word “statistics” should not be taken to imply that every NLME model reduces to a small universal list of moments.
Worked Example: Conditional Moments in a Gaussian Random-Intercept Model
Worked Example: Conditional Moments in a Gaussian Random-Intercept Model
Consider a hypothetical subject with two measurements from the model \(Y_{ij}=\mu+b_i+\epsilon_{ij}\), where \(b_i\sim N(0,\tau^2)\) and independent \(\epsilon_{ij}\sim N(0,\sigma^2)\). At the current parameter values, let \(\mu=4\), \(\tau^2=1\), and \(\sigma^2=1\). The observed measurements are \(y_{i1}=5\) and \(y_{i2}=7\). Find the conditional mean and second moment of \(b_i\), and the expected residual sum of squares needed from this subject.
Assumption check. The random intercept and measurement errors are independent Gaussian variables with the stated variances, and the observations are conditionally independent given \(b_i\). These assumptions give a Gaussian conditional distribution for \(b_i\).
Completing the square in the product of the prior density and the two conditional response densities gives conditional variance \(v_i\) and mean \(m_i\):
Therefore \(b_i\mid y_i;\theta^{(k)}\sim N(4/3,1/3)\). The conditional first and second moments are
For the response likelihood, the residual sum of squares is \(\sum_{j=1}^{2}(y_{ij}-\mu-b_i)^2\). Since \(E[(a-B)^2]=a^2-2aE[B]+E[B^2]\), its conditional expectation is
The subject-specific expected contributions are consequently \(E[b_i\mid y_i]=4/3\), \(E[b_i^2\mid y_i]=19/9\), and an expected residual sum of squares of \(32/9\). These are conditional expectations under the current parameters; this calculation does not update \(\mu\), \(\tau^2\), or \(\sigma^2\).
Why the Conditional Distribution Matters
The posterior weighting follows directly from conditioning. For any integrable function \(h\) of a subject’s random effects,
Substituting the conditional density shows that values of \(b_i\) receive weight in proportion to both their population plausibility and their ability to explain the observed data:
This identity also explains why an unconditional random-effect moment is generally not the E-step answer. The prior expectation \(E[h(B_i)]\) weights effects by the population distribution alone; it ignores the information in \(y_i\).
Worked Example: A Weighted Approximation to a Nonlinear Conditional Expectation
Suppose a hypothetical nonlinear response model has a one-dimensional random effect \(b_i\). At three integration points \(b=-1,0,1\), assume the quadrature contributions to the conditional-density numerator are \(0.10,0.50,0.15\). These contributions already include the integration weights, prior density, and response likelihood. Approximate the conditional expectations of \(b_i\) and \(b_i^2\).
The total contribution is \(0.10+0.50+0.15=0.75\). Dividing each contribution by this total gives normalized weights \(0.10/0.75=2/15\), \(0.50/0.75=2/3\), and \(0.15/0.75=1/5\). They sum to one.
The weighted conditional mean is
The weighted second moment is
These are quadrature approximations, not exact moments unless the chosen integration rule happens to integrate the relevant functions exactly. The nonlinear response affects the relative contributions through \(p(y_i\mid b_i;\theta^{(k)})\); the prior alone would generally produce different weights.
Cross-Moments and Joint Random Effects
For a vector of random effects, it is not enough to calculate each component’s mean and variance separately whenever the complete-data likelihood contains cross-products. For example, a multivariate Gaussian random-effect density involves quadratic terms in the components, including products such as \(b_{i1}b_{i2}\) when the covariance is not diagonal. The relevant conditional expectation is then \(E[b_{i1}b_{i2}\mid y_i]\), not the product of the two conditional means.
Worked Example: A Conditional Cross-Moment
For illustration, suppose numerical integration of a hypothetical subject’s nonlinear likelihood and random-effect density gives the following normalized conditional weights on four possible two-dimensional random-effect values.
| \((b_{i1},b_{i2})\) | Conditional weight | |
|---|---|---|
| \((0,0)\) | 0.10 | |
| \((1,0)\) | 0.20 | |
| \((0,1)\) | 0.30 | |
| \((1,1)\) | 0.40 |
The weights sum to \(1.00\). The conditional means are \(E[b_{i1}\mid y_i]=0.20+0.40=0.60\) and \(E[b_{i2}\mid y_i]=0.30+0.40=0.70\). But the cross-moment is
By contrast, multiplying the means gives \(0.60\times0.70=0.42\), which is not the conditional cross-moment. Using \(0.42\) in place of \(0.40\) would change the expected quadratic terms of a complete-data likelihood that includes \(b_{i1}b_{i2}\).
Numerical Checks, Alternatives, and Limits
In many NLME models, the conditional expectations require numerical integration. The E-step specifies what must be integrated; it does not guarantee that a chosen computation accurately evaluates the integral. Check the calculation at the level of the quantities entering \(Q\), not only at the level of a plotted random-effect estimate.
- Check normalization. Conditional weights should be nonnegative and sum to one after normalization, up to numerical rounding.
- Check refinement stability. Recompute important expected statistics using more quadrature points or a more accurate integration setting. Material changes indicate that the current approximation may be inadequate.
- Check tail coverage. If the likelihood times prior has meaningful mass beyond the numerical integration region, moments—especially second moments—can be biased even when the central region is well represented.
- Check all required terms. Accurate conditional means do not guarantee accurate second moments or cross-moments. Inspect each quantity used in the expected complete-data log-likelihood.
- Check numerical scale. In weighted calculations, very small likelihood contributions can underflow. Stable computation of normalized weights is important when contributions vary greatly.
The exact E-step uses the full conditional distribution and exact conditional expectations. Numerical quadrature approximates those expectations deterministically; simulation-based methods approximate them using draws from a suitable conditional distribution. A posterior mode or empirical Bayes estimate is a different approximation: it substitutes one effect value rather than averaging over conditional uncertainty. It can be useful for other purposes, but it is not the E-step expectation.
The choice among exact calculation, quadrature, and simulation depends on the model’s dimension, shape, and computational demands. No finite set of moments is automatically adequate for every nonlinear model. If the complete-data log-likelihood contains a nonlinear function of \(b_i\), approximating only \(E[b_i]\) and \(E[b_i^2]\) will not in general determine the expected value of that function.
Common Mistake
A common error is to maximize the complete-data log-likelihood after inserting a subject’s conditional mode, then call that calculation an E-step. The mode is one point, whereas the E-step averages over the entire conditional distribution. In general,
where \(\widehat b_i\) is the conditional mode. The two may be close when the conditional distribution is sharply concentrated and the relevant terms vary little over its support, but that is an approximation requiring justification—not the definition of the E-step. A related error is to calculate a nonlinear function of a conditional mean and assume it equals the conditional expectation of that function; generally \(E[h(B_i)\mid y_i]\ne h(E[B_i\mid y_i])\).
Check Your Understanding
Answer each question, then compare with the solutions below.
- In the E-step, which parameter value defines the conditional distribution of the random effects: the candidate value \(\theta\) or the current value \(\theta^{(k)}\)?
- Three normalized conditional weights are \(0.25, 0.50, 0.25\) at random-effect values \(-2,0,2\). Calculate the approximate conditional mean and second moment.
- Why can \(E[b_{i1}b_{i2}\mid y_i]\) not generally be replaced by \(E[b_{i1}\mid y_i]E[b_{i2}\mid y_i]\)?
- State one diagnostic for a quadrature approximation to an E-step expectation and explain what a failure would suggest.
Solutions. 1. The conditional distribution is defined using the current value \(\theta^{(k)}\); \(\theta\) is the candidate argument of \(Q(\theta\mid\theta^{(k)})\). 2. The mean is \((-2)(0.25)+0(0.50)+(2)(0.25)=0\). The second moment is \(4(0.25)+0+4(0.25)=2\). 3. The product of conditional means omits conditional covariance: \(E[b_{i1}b_{i2}\mid y_i]=\operatorname{Cov}(b_{i1},b_{i2}\mid y_i)+E[b_{i1}\mid y_i]E[b_{i2}\mid y_i]\). 4. Increase the quadrature resolution and compare the required expectations. Material changes suggest the original numerical integration did not adequately resolve the conditional distribution or its tails.