The EM View of Latent Random Effects

Why View Random Effects as Latent Data?

In a nonlinear mixed-effects model, each subject has a subject-specific random effect that helps describe that subject’s observations. The random effect is not directly observed; it is integrated out when forming the likelihood for the population parameters. This creates a useful alternative viewpoint: treat the unobserved random effects as latent data, and regard the observations and latent effects together as a complete-data pair.

That viewpoint does not change the statistical model or add measurements. It reorganizes the likelihood calculation. The observed-data likelihood requires integration over random effects, while an expectation-maximization (EM) procedure uses their conditional distribution given the observed data to construct an objective for updating parameters.

Definition: For subject \(i\), let \(y_i\) denote the observed data and \(b_i\) the unobserved random effect. Let \(\theta\) collect the population parameters, including parameters governing the conditional observation model and the random-effect distribution. The joint latent-data density is \(p(y_i,b_i;\theta)\); the observed-data density is its marginal \(p(y_i;\theta)=\int p(y_i,b_i;\theta)\,db_i\). The conditional distribution \(p(b_i\mid y_i;\theta)\) describes the random effect given that subject’s observations at parameter value \(\theta\).

From the NLME Model to the Marginal Likelihood

A standard formulation specifies a conditional observation density \(f(y_i\mid b_i;\theta)\) and a random-effect density \(g(b_i;\theta)\). Their product is the joint density:

p(y_i,b_i;\theta)=f(y_i\mid b_i;\theta)\,g(b_i;\theta).

If subjects are independent under the model, the joint density for all subjects and their random effects is the product of these subject-level densities. The likelihood based only on the observed data is therefore

L(\theta;y)=\prod_{i=1}^{N}p(y_i;\theta) =\prod_{i=1}^{N}\int f(y_i\mid b_i;\theta)\,g(b_i;\theta)\,db_i, \qquad \ell(\theta;y)=\log L(\theta;y).

The integral is essential: the likelihood averages over the model’s possible random effects rather than substituting an estimated effect for each subject. In a nonlinear model, the conditional density can vary sharply with \(b_i\), making these integrals difficult to calculate accurately. Direct marginal-likelihood methods, including numerical integration, address that difficulty by evaluating the integrals.

For EM, instead consider the joint density before integrating out the latent effects. If \(b=(b_1,\ldots,b_N)\), the joint density is \(p(y,b;\theta)\). The complete-data log-likelihood is \(\log p(y,b;\theta)\), viewed as a function of \(\theta\) when both \(y\) and \(b\) are available. In practice \(b\) is unavailable, so EM averages this log density over the conditional distribution of \(b\) given \(y\) at the current parameter value.

The Conditional Random-Effect Distribution and the EM Objective

Bayes’ rule gives the subject-level conditional density:

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

This distribution combines the population distribution of effects with the information about \(b_i\) in subject \(i\)’s observations. It is evaluated at the current parameter value during the expectation step. Write the current estimate as \(\theta^{(t)}\). The EM objective for a candidate parameter \(\theta\) is

Q(\theta\mid\theta^{(t)}) = \operatorname{E}_{b\mid y;\theta^{(t)}} \!\left[\log p(y,b;\theta)\right].

The expectation is over the conditional distribution of all the latent random effects given the observed data under \(\theta^{(t)}\). An EM iteration first evaluates this conditional expectation as a function of the candidate \(\theta\), then chooses a new parameter value that increases or maximizes \(Q(\theta\mid\theta^{(t)})\). The population parameters are updated; the random effects are integrated or averaged over, not treated as newly observed data.

The key justification is a likelihood-ratio identity. A convenient sufficient support condition is that, for the observed data under consideration, the candidate joint density assigns zero probability density wherever the current joint density is zero. In measure-theoretic terms, the candidate joint distribution must be absolutely continuous with respect to the current joint distribution. Also require finite, well-defined likelihoods and the expectations used below. Then

\frac{L(\theta;y)}{L(\theta^{(t)};y)} = \operatorname{E}_{b\mid y;\theta^{(t)}} \left[ \frac{p(y,b;\theta)}{p(y,b;\theta^{(t)})} \right].

To see why, multiply the ratio inside the expectation by the current conditional density and integrate over \(b\). The current joint density cancels, leaving the candidate joint density divided by \(L(\theta^{(t)};y)\). The support condition ensures that this integration does not omit any candidate mass outside the current support:

\int \frac{p(y,b;\theta)}{p(y,b;\theta^{(t)})} \frac{p(y,b;\theta^{(t)})}{L(\theta^{(t)};y)}\,db = \frac{\int p(y,b;\theta)\,db}{L(\theta^{(t)};y)} = \frac{L(\theta;y)}{L(\theta^{(t)};y)}.

Taking logarithms and applying Jensen’s inequality gives a lower bound:

\ell(\theta;y)-\ell(\theta^{(t)};y) \geq \operatorname{E}_{b\mid y;\theta^{(t)}} \left[ \log\frac{p(y,b;\theta)}{p(y,b;\theta^{(t)})} \right] = Q(\theta\mid\theta^{(t)})-Q(\theta^{(t)}\mid\theta^{(t)}).

This is the bridge from the latent-data formulation to EM: a candidate that increases the expected joint log density also increases this lower bound on the log-likelihood change. This derivation requires the support and expectation conditions stated above; the identity is not automatic for arbitrary parameter-dependent supports.

Worked Examples

Worked Example: A Nonlinear Concentration Model and Posterior Effect Weights

Consider one hypothetical subject with one concentration observation. Let the conditional mean be \(C(b)=10\exp[-0.5\exp(b)]\), and suppose the observation has a normal density with standard deviation \(1\). To illustrate marginalization, approximate the random-effect distribution by three quadrature nodes \(b=-0.5,0,0.5\) with weights \(0.25,0.50,0.25\). The observed concentration is \(y=6\). Calculate the approximate marginal density and the posterior weights of the nodes.

At the three nodes, the predicted means are \(C(-0.5)=10\exp(-0.3033)\approx7.384\), \(C(0)=10\exp(-0.5)\approx6.065\), and \(C(0.5)=10\exp(-0.8244)\approx4.385\). For a normal density with standard deviation \(1\), the conditional density at a node is \(\phi(y;C(b),1)=(2\pi)^{-1/2}\exp[-(y-C(b))^2/2]\). Thus the approximate density contributions before multiplying by node weights are about \(0.1531\), \(0.3981\), and \(0.1083\).

The approximate marginal density is the weighted sum:

$$ 0.25(0.1530)+0.50(0.3981)+0.25(0.1083) =0.2644. $$

The corresponding posterior node weights, obtained by normalizing the unrounded weighted contributions, are approximately \(0.1448\), \(0.7528\), and \(0.1024\). They sum to one up to rounding. The observed value gives the node near \(b=0\) the largest posterior weight, but the marginal likelihood still averages over all three nodes. These are quadrature-approximation results, not an exact continuous-distribution calculation.

Worked Example: Conditional Random Effects in a Gaussian Model

For a simple hypothetical calculation, let \(b\sim N(0,1)\) and \(Y\mid b\sim N(b,1)\). An observation is \(y=2\). Find the conditional distribution of \(b\) given \(y\), and compare the marginal density with the conditional density at \(b=2\).

The conditional density for \(b\), up to a constant, is the product of the observation and prior densities: \(p(b\mid y=2)\propto\exp\{-\tfrac12[b^2+(2-b)^2]\}\). Completing the square gives \(b^2+(2-b)^2=2[(b-1)^2+1]\). Therefore the part depending on \(b\) is \(\exp[-(b-1)^2]\), which is the kernel of \(N(1,0.5)\). Hence \(b\mid y=2\sim N(1,0.5)\).

The marginal distribution is \(Y\sim N(0,2)\), so its density at \(2\) is \((4\pi)^{-1/2}\exp(-4/4)\approx0.1038\). By contrast, the conditional observation density evaluated at \(b=2\) is \((2\pi)^{-1/2}\exp[-(2-2)^2/2]\approx0.3989\).

The two numbers answer different questions. The marginal density averages the conditional observation density over the prior distribution of \(b\); the posterior distribution \(N(1,0.5)\) describes plausible effects after observing \(y\). Substituting the posterior mean or a best-fitting effect into the conditional density would not, by itself, produce the marginal likelihood.

Worked Example: Why the Support Condition Matters

Suppose \(y\) has density \(1\) under every parameter value and is independent of \(b\). At the current parameter, let \(b\sim\operatorname{Uniform}(0,1)\); at a candidate parameter, let \(b\sim\operatorname{Uniform}(0,2)\). Compare the marginal likelihood ratio with the expectation of the joint-density ratio under the current conditional distribution.

Because \(y\) has density \(1\) in both cases, both observed-data likelihoods are \(1\); their ratio is \(1\). Under the current parameter, \(b\mid y\) is uniform on \((0,1)\). On that interval, the current joint density is \(1\) and the candidate joint density is \(1/2\). Thus the joint-density ratio is \(1/2\) throughout the current support, and its current conditional expectation is \(1/2\), not \(1\).

The candidate assigns positive probability to \(1<b<2\), where the current density is zero. The expectation under the current conditional distribution cannot include that candidate mass. Therefore the likelihood-ratio identity cannot be used without an appropriate support condition. This example also shows why “integrate the ratio under the current posterior” is not a universally valid shortcut.

Assumptions, Diagnostics, and Alternatives

The latent-data formulation relies on a coherent probability model: the conditional observation density and random-effect density must be normalized, and their product must define the stated joint model. The product across subjects additionally assumes subject-level independence. Any within-subject dependence must be represented in \(f(y_i\mid b_i;\theta)\), rather than ignored. For the likelihood-ratio argument, the candidate and current joint distributions must satisfy the support condition, and the ratios and log expectations must exist as finite quantities.

In computation, a central difficulty is the conditional expectation defining \(Q\). In a nonlinear mixed model, the conditional distribution of \(b\mid y;\theta^{(t)}\) may have no convenient closed form. Numerical integration or simulation may be needed. Useful warning signs include unstable estimates of the conditional expectation when integration settings change, a few quadrature nodes or simulated effects dominating the calculation, and large changes in the result under modest refinements. These signals indicate that the expectation may not be adequately evaluated; they are not evidence that the latent-data identity itself is incorrect.

EM and direct marginal-likelihood maximization use the same observed-data likelihood but organize the computation differently. Direct methods calculate the marginal likelihood, often using numerical integration, and optimize it. EM instead forms a conditional expectation of the joint log density and updates parameters through that objective. EM can make parameter updates easier when the expected joint log density has a convenient form, but its expectation can itself be computationally demanding, especially for nonlinear models. The choice is therefore computational, not a choice between different statistical models.

Common Pitfall

Wrong reasoning: “The random effects are unobserved, so I can replace each one with its current estimate and then maximize the joint density as if those effects were known.” A point estimate does not account for uncertainty in the conditional distribution of \(b\mid y;\theta^{(t)}\); this substitution generally changes the objective and is not the EM expectation step.

Correction: EM averages the joint log density over the full current conditional distribution of the latent effects. Also check support before using the likelihood-ratio identity: if a candidate parameter assigns mass where the current joint density is zero, the expectation under the current conditional distribution can miss that mass.

Key takeaway: The observed-data likelihood integrates out subject-specific random effects. The EM viewpoint retains them as latent variables, conditions on the observed data to describe their uncertainty, and uses the expected joint log density as its parameter-update objective.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

In an NLME model, what distinguishes the joint latent-data density from the observed-data density?
Why does \(p(b_i\mid y_i;\theta^{(t)})\) generally differ from the population random-effect density \(g(b_i;\theta^{(t)})\)?
If two posterior quadrature nodes have normalized weights \(0.2\) and \(0.8\), and a quantity \(r(b)\) equals \(3\) and \(1\) at those nodes, respectively, what is its posterior-weighted expectation?
What support issue invalidates the likelihood-ratio identity in the uniform example?
Name one reason an EM expectation may be difficult to evaluate in a nonlinear mixed-effects model.

Solutions. (1) The joint density includes both \(y_i\) and \(b_i\); the observed-data density integrates \(b_i\) out. (2) The observations update the relative plausibility of effects, so the conditional distribution combines the observation density with the population distribution. (3) \(0.2(3)+0.8(1)=1.4\). (4) The candidate distribution has positive mass outside the support of the current joint distribution, so the current conditional expectation omits candidate mass. (5) The conditional random-effect distribution or its expected log density may lack a closed form and require accurate numerical integration or simulation.