Why Exact EM Is Difficult in Nonlinear Mixed Models

Why the E-Step Becomes Difficult

In a nonlinear mixed-effects model, each individual’s observations share latent random effects. The E-step of EM asks for expectations of complete-data quantities conditional on the observed data and the current parameter values. The difficulty is not that these expectations are undefined: under ordinary model conditions they are well-defined integrals. The difficulty is that the conditional distribution inside those integrals is often not a standard distribution, and the integrals generally have no closed-form solution.

The contrast with a linear-Gaussian model is important. There, the conditional distribution of the random effects is Gaussian, so moments such as \(E[b_i\mid y_i]\) and \(E[b_i^2\mid y_i]\) can be calculated directly. With a nonlinear individual prediction, the likelihood as a function of the random effects usually has a non-quadratic log density. Multiplying it by a Gaussian random-effects density does not, in general, yield another Gaussian distribution.

Definition: An exact E-step calculates the conditional expectations of the complete-data quantities required by the M-step, integrating over the conditional distribution of the latent random effects given the observations and the current parameter values. “Exact” refers to evaluating these expectations without numerical or stochastic approximation; it does not mean that the resulting parameter estimates are exact for the population.

For individual \(i\), let \(y_i\) denote the observation vector, \(b_i\) the random-effects vector, and \(\theta\) the population parameters. Suppose the model specifies a conditional density \(p(y_i\mid b_i,\theta)\) and a random-effects density \(p(b_i\mid\theta)\). The conditional density needed by the E-step is

p(b_i\mid y_i,\theta)= \frac{p(y_i\mid b_i,\theta)p(b_i\mid\theta)} {\int p(y_i\mid u,\theta)p(u\mid\theta)\,du}.

For a complete-data quantity \(S_i(b_i,y_i;\theta)\), its conditional expectation is therefore

E[S_i\mid y_i,\theta]= \frac{\int S_i(b_i,y_i;\theta)\,p(y_i\mid b_i,\theta)p(b_i\mid\theta)\,db_i} {\int p(y_i\mid b_i,\theta)p(b_i\mid\theta)\,db_i}.

The denominator is the individual’s marginal likelihood contribution. Even when it can be evaluated numerically, the numerator may require a separate integral for each quantity needed in the M-step. The earlier material on the E-step as conditional sufficient-statistic calculation explains why these expectations, rather than just a single estimate of \(b_i\), are needed.

The Linear-Gaussian Case as a Contrast

First consider a random-intercept model with known observation variance. For one individual, suppose \(y_{ij}=b_i+\epsilon_{ij}\), where \(b_i\sim N(0,\omega^2)\) and the residuals are conditionally independent with \(\epsilon_{ij}\sim N(0,\sigma^2)\). The observations are conditionally independent given \(b_i\), and the random effect is independent of the residuals. The log conditional density of \(b_i\), up to a constant, is

-\frac{b_i^2}{2\omega^2} -\sum_{j=1}^{n_i}\frac{(y_{ij}-b_i)^2}{2\sigma^2}.

This is a quadratic function of \(b_i\). Completing the square gives a Gaussian conditional distribution, from which the required moments follow. The quadratic form is the feature that makes the calculation tractable, not simply the fact that the prior random effect is Gaussian.

Worked Example: A Closed-Form Random-Intercept E-Step

Suppose one hypothetical individual has three observations \(1.2,\ 1.4,\ 1.6\), with \(b\sim N(0,1)\) and residual variance \(\sigma^2=0.25\). Assume conditional independence and treat the variances as known for this E-step.

The conditional precision is the prior precision plus the observation precision contributed by each measurement: \(1+3/0.25=13\). Thus the conditional variance is \(1/13=0.07692\). The conditional mean is the precision-weighted mean:

E[b\mid y]= \frac{\sum_j y_j/0.25}{1+3/0.25} =\frac{(1.2+1.4+1.6)/0.25}{13} =\frac{16.8}{13} =1.2923.

The second moment needed for a variance update is the conditional variance plus the square of the conditional mean: \(E[b^2\mid y]=0.07692+(1.2923)^2\approx1.747\). These moments are available from the Gaussian conditional distribution. This example establishes a useful baseline: an E-step can be analytically straightforward when the complete-data model is conjugate.

How Nonlinear Predictions Change the Conditional Density

Now suppose observations are normally distributed around a nonlinear individual prediction \(f(t_{ij};\theta,b_i)\), with residual variance \(\sigma^2\), and \(b_i\sim N(0,\Omega)\). For one individual, the log conditional density of \(b_i\), apart from constants, is

-\frac{1}{2}b_i^{T}\Omega^{-1}b_i -\frac{1}{2\sigma^2}\sum_j \left[y_{ij}-f(t_{ij};\theta,b_i)\right]^2.

If \(f\) is nonlinear in \(b_i\), the squared residual term generally is not quadratic in \(b_i\). Consequently, the log conditional density is not generally quadratic, and the conditional density is not generally Gaussian. The E-step may require expectations of nonlinear functions of \(b_i\), not merely a conditional mean. In particular, expanding the squared residual can produce terms involving \(f\) and \(f^2\), whose conditional expectations must also be calculated.

Worked Example: An Exponential Individual Prediction

Consider one hypothetical observation \(y=2\) from \(y=\exp(b)+\epsilon\), where \(b\sim N(0,1)\) and \(\epsilon\sim N(0,0.25)\). The fixed-effect part of the prediction is set to zero for simplicity. Assume the observation is conditionally normal and independent of the random effect before conditioning on the data.

The conditional density kernel for \(b\) is proportional to \(\exp\{-b^2/2-2(2-\exp(b))^2\}\). Its log kernel is \(-b^2/2-2(2-\exp(b))^2\). Differentiating twice with respect to \(b\) gives

\frac{d^2}{db^2}\log p(b\mid y) =-1+8\exp(b)-8\exp(2b).

For a Gaussian density, the second derivative of its log density is constant in \(b\). Here it equals \(-1\) at \(b=0\), but equals \(-17\) at \(b=\log 2\). Therefore the conditional density is not Gaussian. The complete-data residual term is \((2-\exp(b))^2=4-4\exp(b)+\exp(2b)\), so an E-step involving this term requires conditional expectations of \(\exp(b)\) and \(\exp(2b)\). The Gaussian prior alone does not supply those conditional expectations: they must be computed under the non-Gaussian posterior density.

This example illustrates the general obstacle. The normalization integral for \(p(b_i\mid y_i,\theta)\) and the expectation integrals for the complete-data quantities contain exponentials of nonlinear functions. In common pharmacometric predictions—such as exponential decay or compartmental concentration functions—these integrals need not reduce to a familiar distribution’s moments. Special parameter settings or model forms can still yield closed forms, so “nonlinear” means that a closed form is generally unavailable, not that it is impossible in every case.

A Numerical Approximation to a Nonlinear E-Step Moment

A deterministic integration rule can approximate an expectation by evaluating its numerator and denominator integrands at a finite set of points. The next example uses composite Simpson quadrature. It demonstrates the mechanics, not a recommended final accuracy: the integration range is deliberately limited so that the calculation is transparent. In practice, both range and resolution must be checked.

Worked Example: Approximating a Second Moment by Quadrature

Suppose a hypothetical observation is \(y=1\) from the nonlinear model \(y=b^2+\epsilon\), with \(b\sim N(0,1)\) and \(\epsilon\sim N(0,1)\). The assumptions are conditional normality and independence of \(b\) and \(\epsilon\). For estimating the random-effect variance, the E-step needs \(E[b^2\mid y]\). Ignoring constants common to all \(b\), the posterior kernel is \(k(b)=\exp\{-b^2/2-(1-b^2)^2/2\}\).

Apply composite Simpson quadrature on \([-2,2]\) with step \(h=0.5\), using nodes \(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2\). The corresponding kernel values, rounded to four significant figures, are \(0.001503,\ 0.1486,\ 0.6065,\ 0.6661,\ 0.6065,\ 0.6661,\ 0.6065,\ 0.1486,\ 0.001503\). For example, at \(b=0.5\), \(k(0.5)=\exp\{-0.5(0.25)-0.5(0.75)^2\}=\exp(-0.40625)\approx0.6661\).

The Simpson weights are \(1,4,2,4,2,4,2,4,1\). Applying them to the kernel values gives a weighted sum of approximately \(10.160\), so the common factor \(h/3\) gives a denominator integral of approximately \(1.693\). For the numerator, multiply each kernel value by \(b^2\) before applying the same weights. The weighted sum is approximately \(6.445\), giving a numerator integral of approximately \(1.074\). Therefore

E[b^2\mid y]\ \approx\ \frac{\int_{-2}^{2}b^2k(b)\,db}{\int_{-2}^{2}k(b)\,db} \approx\frac{1.074}{1.693} =0.6344.

The ratio is an approximation to the truncated integrals, and hence to the desired full-line expectation; it is not an exact E-step result. The omitted tails and the quadrature discretization both contribute error. The calculation shows how even a single required moment becomes an integration problem when the conditional density is non-Gaussian.

Checking Approximations and Choosing an Approach

When a closed-form conditional expectation is unavailable, there are several alternatives. Deterministic quadrature evaluates integrals at selected nodes; importance sampling estimates them using weighted draws; Laplace methods approximate the conditional density locally near its mode; and Monte Carlo methods approximate expectations with simulated random effects. These approaches differ in cost, stability, and approximation error. Earlier material on numerical integration in NLME models treats quadrature, importance sampling, and integration diagnostics in detail. The key point here is that these methods approximate the E-step expectation; they do not turn a non-Gaussian conditional distribution into a Gaussian one.

A numerical E-step should be checked rather than trusted solely because an algorithm returns a value. Useful checks include increasing the number of quadrature nodes, expanding the integration range, comparing a deterministic calculation with a separate integration strategy, and confirming that the calculated expectations are stable at the precision needed by the M-step. The small-range calculation above is especially sensitive to tail omission; agreement after refining nodes over the same narrow range would not by itself establish that the tails are negligible.

A common mistake is to replace \(E[b_i\mid y_i]\) by a single best-fitting random effect and treat that plug-in value as the E-step. In general, \(E[S(b_i)\mid y_i]\) is not \(S(E[b_i\mid y_i])\). For instance, even when a conditional mean is known, \(E[b_i^2\mid y_i]\) also depends on the conditional variance. For nonlinear functions, the discrepancy can be larger: \(E[\exp(b_i)\mid y_i]\) cannot generally be obtained by exponentiating \(E[b_i\mid y_i]\). The E-step requires the conditional expectation of each complete-data quantity specified by the model, not merely a point estimate of the latent effect.

Key takeaway: Exact EM is analytically difficult in nonlinear mixed models because the E-step integrates complete-data quantities against a conditional random-effects distribution that is usually non-Gaussian. Closed forms remain possible in special cases, but otherwise expectations must be approximated and their numerical accuracy assessed.

Check Your Understanding

Answer each question, then compare with the brief solutions.

Why does a linear-Gaussian random-intercept model permit a closed-form conditional mean for the random effect?
In the exponential-prediction example, why is the conditional distribution not Gaussian?
For a variance update, why is \(E[b^2\mid y]\) not generally equal to \(\{E[b\mid y]\}^2\)?
Name two checks that can reveal numerical integration error in an approximate E-step.
Does the nonlinearity of a prediction prove that no closed-form E-step can exist? Explain.

Solutions. (1) The log conditional density is quadratic in the random effect, so completing the square gives a Gaussian conditional distribution with calculable moments. (2) Its log density has a second derivative that varies with \(b\), whereas a Gaussian log density has constant second derivative. (3) The second moment equals the conditional variance plus the square of the conditional mean, so it includes uncertainty beyond the point estimate. (4) Refine the quadrature nodes, expand the integration range, or compare with another integration method. (5) No. Nonlinearity makes a closed form generally unavailable, but special forms or parameter settings may still permit one.