Conditional Random-Effect Simulation in SAEM

Why Condition on the Observed Data?

In a mixed-effects model, a subject’s random effects are latent: they are not observed directly, but they help explain that subject’s measurements. During SAEM, the E-step requires information about these latent effects conditional on the data actually observed and on the current population-parameter values. Drawing random effects from their population distribution alone does not provide that information: it ignores how well each possible effect explains the subject’s measurements.

The target distribution is the conditional distribution of a subject’s random effects given that subject’s observations. In a nonlinear model it is generally not a familiar distribution and may not be available in closed form. SAEM can nevertheless use draws from this conditional distribution to simulate the complete-data statistics needed for its stochastic E-step.

Definition: Let $y_i$ denote all observations for subject $i$, $\eta_i$ that subject’s random-effect vector, and $\theta$ the population and residual parameters. At iteration $k$, the conditional random-effect distribution is $p(\eta_i\mid y_i,\theta^{(k)})$. Under the model’s usual factorization, its density is proportional to the subject-specific likelihood times the random-effect density: $$ p(\eta_i\mid y_i,\theta^{(k)}) = \frac{p(y_i\mid \eta_i,\theta^{(k)})p(\eta_i\mid\theta^{(k)})} {\int p(y_i\mid u,\theta^{(k)})p(u\mid\theta^{(k)})\,du}. $$ The denominator normalizes the density over the random-effect space.

Here $p(\eta_i\mid\theta^{(k)})$ is the population distribution of random effects, while $p(y_i\mid\eta_i,\theta^{(k)})$ is the likelihood for that subject conditional on the random effects. The conditioning on $\theta^{(k)}$ matters: as population or residual parameters change, the conditional distribution used for the next SAEM E-step also changes.

From the Conditional Distribution to the SAEM E-Step

Let $T(y,\eta)$ be the vector of complete-data sufficient statistics used by the model’s M-step. In exact EM, the E-step evaluates the conditional expectation of these statistics given all observations and the current parameters. In SAEM, the expectation is approximated using simulated random effects from the conditional distributions. For example, if a draw $\eta_i^{(k)}$ is obtained for each subject, the simulated statistic is $\sum_i T_i(y_i,\eta_i^{(k)})$. The stochastic approximation then combines this simulated statistic with the running statistic, rather than treating one draw as an exact conditional expectation.

The subject-level conditioning follows from the model factorization: given the population parameters and the random effects, subjects’ observation vectors are independent, and their random effects are independent draws from the population distribution. Consequently, the joint conditional distribution factors across subjects:

p(\eta_1,\ldots,\eta_N\mid y_1,\ldots,y_N,\theta^{(k)}) = \prod_{i=1}^{N}p(\eta_i\mid y_i,\theta^{(k)}).

This factorization permits subject-by-subject conditional simulation. It does not mean that the components within a subject’s random-effect vector are independent: if the population covariance matrix has off-diagonal elements, the conditional components will generally remain dependent.

Worked Example: A Gaussian Conditional Distribution

Suppose a subject has one measurement $y=1.20$ from the model $y\mid b\sim N(b,0.25)$, where $b\sim N(0,1)$. Treat the population and residual parameters as fixed at their current values. Find the conditional distribution of $b$ given $y$.

Ignoring factors that do not depend on $b$, the likelihood times prior density is proportional to $\exp[-(1.20-b)^2/0.50-b^2/2]$. Expanding the exponent gives

-\frac{(1.20-b)^2}{0.50}-\frac{b^2}{2} =-2.5b^2+4.8b-2.88 =-\frac{1}{2}\left(5b^2-9.60b+5.76\right).

The constant term inside the parentheses is $5.76$. Completing the square, $5b^2-9.60b+5.76=5(b-0.96)^2+1.152$. The part depending on $b$ is therefore $-\frac{5}{2}(b-0.96)^2$, the exponent of a normal density with mean $0.96$ and variance $1/5=0.20$. Thus

b\mid y=1.20\sim N(0.96,0.20).

This result also follows from Gaussian precision addition: the prior precision is $1$, the observation precision is $1/0.25=4$, and the posterior precision is $5$. If a standard-normal draw $z=-0.50$ is used to generate a conditional draw, then $b=0.96+\sqrt{0.20}(-0.50)=0.736$ (to three significant figures). This is one simulated value, not the conditional mean. The conditional mean remains $0.96$, and the conditional expectation of $b^2$ is $0.20+0.96^2=1.122$.

Nonlinear Likelihoods and Conditional Weights

For a nonlinear observation model, Bayes’ formula still defines the conditional distribution, but completing the square will not generally yield a normal distribution. For instance, with a nonlinear prediction $f(t_{ij},\eta_i)$, conditionally independent Gaussian residuals with variance $\sigma^2$, and random effects $\eta_i\sim N(0,\Omega)$, the conditional density is proportional to

\exp\left\{ -\frac{1}{2\sigma^2}\sum_j \bigl(y_{ij}-f(t_{ij},\eta_i)\bigr)^2 -\frac{1}{2}\eta_i^{\mathsf T}\Omega^{-1}\eta_i \right\}.

The first term rewards random effects that fit the subject’s observations; the second reflects their plausibility under the population distribution. Both terms matter. A random effect may fit the observations well but still have little conditional support if it is far into the population distribution’s tail.

Worked Example: Relative Conditional Weights in an Exponential Model

Consider one hypothetical measurement $y=0.80$ with prediction $f(b)=\exp(-b)$, residual distribution $y\mid b\sim N(\exp(-b),0.04)$, and prior $b\sim N(0,0.25)$. The conditional density, apart from a common normalizing constant, is

p(b\mid y)\ \propto\ \exp\left[-\frac{(0.80-\exp(-b))^2}{0.08}\right] \exp\left[-\frac{b^2}{0.50}\right].

At $b=-0.50$, the prediction is $\exp(0.50)=1.64872$. The likelihood factor and prior factor are

\exp\left[-\frac{(0.80-1.64872)^2}{0.08}\right] \approx 0.0001229,\qquad \exp\left[-\frac{(-0.50)^2}{0.50}\right] =\exp(-0.50)\approx 0.6065.

Their product is an unnormalized conditional weight of approximately $0.0001229066\times0.6065307=0.00007455$. At $b=0$, the likelihood factor is $\exp[-(0.80-1)^2/0.08]=\exp(-0.50)\approx0.6065$, and the prior factor is $1$, so the unnormalized weight there is approximately $0.6065$. The much smaller weight at $-0.50$ reflects its poor fit to this observation, despite its nonzero prior density.

These pointwise weights illustrate the two factors in the conditional density; they are not posterior probabilities at the grid points. A finite grid can help visualize or approximate a distribution, but a draw intended for SAEM must represent the full conditional distribution, not simply select a point with a large unnormalized density.

Using Draws to Build the Stochastic E-Step

A conditional draw supplies a simulated value of the latent part of the complete data. Substituting that draw into the complete-data statistics creates a random approximation to the conditional expectation required by exact EM. Because it is random, one draw can differ appreciably from the expectation; SAEM handles this by stochastic approximation across iterations. The key requirement for this tutorial is that the draws target $p(\eta_i\mid y_i,\theta^{(k)})$, not the prior distribution.

Worked Example: Two Subject Draws and a Sufficient Statistic

Suppose two hypothetical subjects each have one observation in a Gaussian random-intercept model $y_i\mid b_i\sim N(b_i,1)$, with $b_i\sim N(0,1)$. The current parameters are fixed, and the observed values are $y_1=1$ and $y_2=-1$. For each subject, the conditional variance is $1/(1+1)=0.50$, and the conditional mean is $y_i/2$. Hence $b_1\mid y_1\sim N(0.50,0.50)$ and $b_2\mid y_2\sim N(-0.50,0.50)$.

For illustration, take independent standard-normal simulation values $z_1=0.40$ and $z_2=-1.00$. The resulting conditional draws are

b_1^{(k)}=0.50+\sqrt{0.50}(0.40)=0.7828,\qquad b_2^{(k)}=-0.50+\sqrt{0.50}(-1.00)=-1.2071.

If a complete-data statistic is $T=\sum_i b_i^2$, this simulation gives $T^{(k)}=(0.7828)^2+(-1.2071)^2\approx0.6128+1.4571=2.0699$. The exact conditional expectation is different: $E(T\mid y)=\sum_i\{\operatorname{Var}(b_i\mid y_i)+[E(b_i\mid y_i)]^2\} =(0.50+0.25)+(0.50+0.25)=1.50$.

The difference between $2.0699$ and $1.50$ is ordinary simulation variation, not an arithmetic contradiction. The simulated statistic is an input to the SAEM stochastic approximation; it is not asserted to equal the exact conditional expectation on every iteration.

Assumptions, Diagnostics, and Alternatives

The conditional-density expressions above rely on the stated probability model. In the Gaussian examples, the residuals are conditionally normal with the specified variances, and the random effects have the stated normal prior distributions. The subject-level factorization additionally assumes independent subjects under the model. In an application, the conditional likelihood must include all relevant observations for that subject and use the current values of the model parameters. Omitting observations, using stale parameter values, or using the wrong variance changes the target distribution.

A useful diagnostic question is whether the proposed simulation mechanism produces values in regions with appreciable conditional density. For a nonlinear model, compare the likelihood and prior contributions over plausible random-effect values, and check that simulated values are not systematically confined to implausible extremes or to a narrow region that misses relevant modes. In multiple dimensions, inspect each component and their joint behavior; marginal summaries alone can hide dependence. Repeated simulations at fixed parameters can also reveal whether empirical summaries are stable enough for the intended stochastic update.

Conditional simulation is not the only way to approximate the E-step. Numerical integration, including quadrature, can approximate conditional expectations directly when the random-effect dimension is manageable and the integrand is adequately resolved. Simulation can be more practical when the integral is difficult or the dimension is larger, but it introduces Monte Carlo variability and depends on obtaining representative conditional draws. A proposal-based simulation method is one possible way to generate such draws when direct sampling is unavailable; the central requirement remains that the target distribution is the conditional distribution, not an arbitrary convenient approximation.

Approach	Strength	Limitation
Direct conditional sampling	Produces draws from the target when a tractable conditional form is available.	Often unavailable for nonlinear mixed-effects models.
Conditional numerical integration	Can approximate conditional expectations without simulation noise when the integral is well resolved.	Computational cost and integration error can grow with dimension or difficult geometry.
Simulation-based approximation	Can represent complicated conditional distributions and generate complete-data statistics.	Finite simulation creates variability; poor exploration gives unrepresentative statistics.

A common mistake is to draw each subject’s random effect from its population distribution and call that the SAEM E-step. That draw ignores the subject’s data: in the exponential example, the prior alone would not penalize $b=-0.50$ for predicting $1.64872$ when the observation is $0.80$. The correction is to target the likelihood-times-prior conditional distribution. A second mistake is to interpret a single simulated statistic as the exact conditional expectation; the two-subject example shows why that need not hold.

Key takeaway: SAEM’s simulated E-step requires random-effect draws conditional on each subject’s observations and the current parameters. Those draws turn latent complete-data statistics into simulated quantities that can be incorporated through stochastic approximation.

Check Your Understanding

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

In one sentence, distinguish the random-effect prior from the conditional random-effect distribution used in the SAEM E-step.
For the Gaussian example with $y=1.20$, state the conditional mean and variance of $b$, and calculate $E(b^2\mid y)$.
In the exponential example, why is the unnormalized weight at $b=-0.50$ not itself a posterior probability?
For the two-subject example, explain why the simulated statistic $2.0699$ can differ from the exact conditional expectation $1.50$.

Solutions. (1) The prior describes random effects before conditioning on the subject’s observations; the E-step uses the posterior conditional distribution, proportional to the subject likelihood times that prior. (2) The conditional distribution is $N(0.96,0.20)$, so $E(b^2\mid y)=0.20+0.96^2=1.1216$. (3) It is an unnormalized density value at one point; a posterior probability requires a normalized probability for an event or region, not a pointwise density. (4) The simulated statistic is based on one random draw per subject, so it varies around the conditional expectation rather than equaling it on every iteration.