Tutorials › Nonlinear Mixed-Effects Models (NLME) › How Many Latent-Effect Draws Per SAEM Iteration

The EM and SAEM Algorithms · Tutorial 212 of 797

How Many Latent-Effect Draws Per SAEM Iteration

Compare one and multiple conditional random-effect draws per SAEM iteration, accounting for Monte Carlo noise, computational cost, and dependence between draws.

Advanced 10 min read

What You'll Learn

  • Define a batch of conditional random-effect draws and its role in a SAEM update.
  • Derive how batch size and autocorrelation affect the variability of a sufficient-statistic estimate.
  • Calculate the effective information in correlated draws.
  • Compare single-draw and batch-average strategies without changing the SAEM gain schedule.
  • Recognize when more draws add little precision or introduce a misleading sense of stability.

Why the Number of Draws Is a Design Choice

In SAEM, a latent random-effect state is simulated conditionally on the observed data and the current population-parameter estimate. Its sufficient statistic contributes to a stochastic approximation, which is then used in the parameter update. A practical choice is how many latent-effect states to use for that contribution at each iteration.

One state per iteration is computationally economical and can be adequate when many SAEM iterations are feasible and the conditional simulation mixes reasonably. Averaging several states can reduce the simulation noise in an iteration, but costs more and may yield much less precision than expected if the states are highly dependent. The relevant quantity is not simply the number of draws; it is the amount of effectively independent information they provide.

Definition: Let \(y\) denote the observed data, \(\theta_{t-1}\) the population parameters entering iteration \(t\), and \(b\) the latent random-effect state. Write \(s(y,b)\) for the sufficient-statistic contribution used by the SAEM update. With \(K\) conditional draws \(b_t^{(1)},\ldots,b_t^{(K)}\), the batch statistic is $$ \bar{s}_t=\frac{1}{K}\sum_{k=1}^{K}s(y,b_t^{(k)}). $$ A batch-average SAEM update is \(S_t=S_{t-1}+\gamma_t(\bar{s}_t-S_{t-1})\), where \(S_{t-1}\) is the running statistic and \(\gamma_t\) is the gain for outer iteration \(t\).

A “draw” here means a latent-effect state contributing to the statistic. Depending on the simulation scheme, that state may include random effects for all subjects or a specified block of them. In a Markov chain implementation, successive states are generally dependent. Also, averaging \(K\) states within one outer iteration is not the same operation as applying \(K\) separate SAEM updates: the batch-average formula applies the gain once.

What Batch Size Does to Monte Carlo Variability

First consider a scalar component \(G=s(y,b)\) at a fixed data set and parameter value. Suppose the simulated states are stationary draws from the desired conditional distribution, so each \(G_k=s(y,b_t^{(k)})\) has conditional mean \(m\) and variance \(\sigma_G^2\). If the \(K\) draws are independent, the variance of their average is \(\sigma_G^2/K\). Thus, increasing \(K\) by a factor of four halves the standard deviation of the batch average, not quarters it.

For dependent draws from a stationary Markov chain, define \(\rho_h\) as the lag-\(h\) autocorrelation of \(G_k\). The variance of the average is then:

$$ \operatorname{Var}(\bar{s}_t)=\frac{\sigma_G^2}{K} \left[1+2\sum_{h=1}^{K-1}\left(1-\frac{h}{K}\right)\rho_h\right]. $$

The bracketed term is the variance inflation relative to \(K\) independent draws. Positive autocorrelation makes that term greater than one, so the variance is reduced less than the independent-draw calculation predicts. For a particular statistic and chain, an approximate effective sample size is \(K\) divided by this inflation factor. It is a way to describe the information in the batch, not a new count of actual simulations.

This variance calculation assumes stationarity for the batch and finite variance for the statistic. A chain that has not reached the intended conditional distribution can also create bias; averaging more states does not, by itself, remove such bias. The calculation is most directly applicable to a fixed conditional target. During SAEM the parameters change, so the target may change as well, and the relationship between within-batch and across-iteration dependence matters.

Worked Example: Independent Draws and the Expected Noise Reduction

Suppose a component of the conditional statistic has variance \(9\) units squared under the current conditional distribution. Compare one draw with a batch of four independent draws.

For one draw, the variance of the statistic contribution is \(9\) units squared, and its standard deviation is \(\sqrt{9}=3\) units.

For four independent draws, the variance of their average is: $$ \operatorname{Var}(\bar{s}_t)=\frac{9}{4}=2.25\text{ units}^2, \qquad \operatorname{SD}(\bar{s}_t)=\sqrt{2.25}=1.5\text{ units}. $$

The four-draw average has half the standard deviation of a single draw. It requires four times as many conditional simulations for that iteration. This is the idealized reduction: it depends on independence and correct sampling from the conditional distribution.

Worked Example: How Autocorrelation Changes the Benefit

Now suppose the same statistic variance is \(9\) units squared, but four successive draws have an approximate first-order autocorrelation pattern \(\rho_h=0.6^h\). Use the finite-batch variance formula to assess the four-draw average.

For \(K=4\), the variance inflation factor is: $$ 1+2\left[\left(1-\frac{1}{4}\right)(0.6) +\left(1-\frac{2}{4}\right)(0.36) +\left(1-\frac{3}{4}\right)(0.216)\right] =1+2(0.450+0.180+0.054)=2.368. $$

Therefore: $$ \operatorname{Var}(\bar{s}_t)=\frac{9}{4}(2.368)=5.328\text{ units}^2, \qquad \operatorname{SD}(\bar{s}_t)=\sqrt{5.328}\approx2.308\text{ units}. $$

The four correlated draws still give a lower variance than one draw under the stationary assumptions, but the standard deviation is about \(2.308\) units rather than the \(1.5\) units predicted by independence. The effective sample size is approximately \(4/2.368=1.689\) independent draws. Counting four simulated states as four independent contributions would overstate the precision.

The SAEM Update and the Gain

Batch size and gain schedule are distinct choices. The batch statistic summarizes the conditional draws within iteration \(t\); the gain \(\gamma_t\) determines how strongly that iteration changes the running statistic. If the intended update uses one gain application per outer iteration, averaging more draws should ordinarily change \(\bar{s}_t\), not silently turn one iteration into \(K\) gain applications.

Worked Example: Applying the Gain Once to a Batch

Suppose the current running statistic is \(S_{t-1}=10\), the gain is \(\gamma_t=0.2\), and three hypothetical conditional draws produce statistic values \(8\), \(12\), and \(11\). Calculate the batch-average update.

The batch average is: $$ \bar{s}_t=\frac{8+12+11}{3}=\frac{31}{3}\approx10.333. $$

Apply the gain once: $$ S_t=10+0.2(10.333-10)=10+0.0667\approx10.067. $$

For comparison, using only the first draw for this iteration would give \(10+0.2(8-10)=9.6\). That difference illustrates the random variation a batch can smooth; it does not mean the batch update is automatically closer to the conditional expectation in every realized iteration. Applying the gain separately to all three values would be a different update rule and should not be confused with the batch-average calculation.

A Short Derivation of the Dependence Adjustment

The variance formula follows by expanding the variance of a sum. Under stationarity, every \(G_k\) has variance \(\sigma_G^2\), and the covariance between \(G_i\) and \(G_j\) is \(\sigma_G^2\rho_{|i-j|}\). Consequently:

$$ \begin{aligned} \operatorname{Var}(\bar{G}) &=\frac{1}{K^2}\operatorname{Var}\left(\sum_{k=1}^{K}G_k\right)\\ &=\frac{1}{K^2}\left[K\sigma_G^2+ 2\sigma_G^2\sum_{h=1}^{K-1}(K-h)\rho_h\right]\\ &=\frac{\sigma_G^2}{K} \left[1+2\sum_{h=1}^{K-1}\left(1-\frac{h}{K}\right)\rho_h\right]. \end{aligned} $$

In the second line, there are \(K\) variance terms and, for each lag \(h\), \(K-h\) pairs separated by that lag. The factor of two accounts for both orders of each pair in the variance expansion. This derivation shows why autocorrelation in the statistic that enters the update is the important dependence to examine. A chain can appear to move in one component while the sufficient statistic of interest remains strongly autocorrelated.

Choosing a Practical Number of Draws

There is no universally optimal \(K\). The choice depends on the cost of conditional simulation, the dependence of the simulated statistics, the gain schedule, and the number of useful outer iterations available. The main options within this choice are:

StrategyPotential advantageMain caution
One contributing draw per iterationLow simulation cost per iteration; many outer updates can be made for a fixed simulation budget.Each statistic update can be noisy, and poor chain mixing can make consecutive updates dependent.
Average several draws per iterationCan reduce the noise in the statistic supplied to that iteration.Costs more; positive chain autocorrelation can substantially reduce the benefit.
Run extra chain transitions but use only a later stateCan allow the state to move farther from its starting value before contributing.This is not the same as averaging the transitions; extra computation does not guarantee adequate mixing.

For the first comparison, a single draw is not necessarily a fresh independent draw from the conditional distribution. In a Markov chain implementation it is usually the current state after a transition, and its dependence on previous states matters. Conversely, a batch of several draws need not be independent to be useful, but its expected precision should be judged using the autocorrelation of the statistic, not just \(K\).

A reasonable implementation decision is therefore empirical and statistic-specific. Examine traces and autocorrelation for relevant components of \(s(y,b)\), not only for the raw random effects. If the batch average is still highly variable, consider whether additional draws meaningfully increase effective information, whether the chain is moving adequately, and whether the conditional target is being sampled as intended. Assess computational cost per useful update as well as time per simulation.

Diagnostics, Failure Modes, and Limitations

More draws can reduce Monte Carlo variability only under appropriate sampling conditions. If a chain has not explored the relevant conditional distribution, a larger batch of near-duplicate states may give a stable-looking average that is nevertheless unrepresentative. When conditional simulations are initialized or updated in a way that creates transient behavior, increasing \(K\) does not automatically remove the resulting bias. The assumptions behind the variance calculation—stationarity and finite variance—must be credible before using it as a precision guide.

  • Check the statistic that is actually averaged. Autocorrelation in a component of \(b\) may differ from autocorrelation in a nonlinear function of \(b\) used by the sufficient statistic.
  • Compare nominal and effective batch size. If positive dependence is strong, a large \(K\) can still carry information similar to a much smaller independent sample.
  • Distinguish a batch average from repeated SA updates. Applying the gain once to an average and applying it repeatedly to individual draws generally produce different updates.
  • Keep the simulation cost in view. A lower-noise update is useful only if its added cost is worthwhile relative to using that computation for further outer iterations or improving the simulation scheme.

These calculations are guides, not a guarantee that a particular batch size will produce a more accurate final parameter estimate. SAEM updates depend on a changing parameter value and a sequence of stochastic statistics; a fixed-target variance calculation isolates only one part of that process. Convergence still relies on suitable conditional simulation and gain behavior. Batch size does not replace those requirements.

Common pitfall: “I used \(K\) draws, so the Monte Carlo standard error must be \(1/\sqrt{K}\) times the one-draw standard error.” That reduction follows only for independent draws with the stated variance. For dependent draws, the finite-batch autocorrelation factor inflates the variance; for draws from a poorly explored or nonstationary chain, even that variance calculation may not describe the error.

Check Your Understanding

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

  1. If independent draws have statistic variance \(16\) units squared, what are the variance and standard deviation of the average of four draws?
  2. In the four-draw example with variance inflation factor \(2.368\), why is the effective sample size less than four?
  3. With \(S_{t-1}=5\), \(\gamma_t=0.1\), and batch statistic values \(4, 6, 8\), calculate the batch-average SAEM update.
  4. Does increasing the number of draws necessarily remove bias caused by a chain that has not reached its intended conditional distribution? Explain.

Solutions. (1) The variance is \(16/4=4\) units squared, and the standard deviation is \(\sqrt{4}=2\) units. (2) Positive autocorrelation makes the variance larger than the independent-draw variance, so four dependent states contain less information than four independent states; \(4/2.368\approx1.689\). (3) The batch average is \((4+6+8)/3=6\), so \(S_t=5+0.1(6-5)=5.1\). The gain is applied once. (4) No. If the chain has not explored the target distribution, the batch may remain unrepresentative; averaging more such draws does not by itself eliminate that bias.