What Reproducibility Means for SAEM
SAEM uses simulation to approximate quantities needed for parameter estimation. Consequently, two runs of the same model can follow different latent-effect draws and produce slightly different parameter traces. Reproducibility requires more than saving the final estimates: another analyst should be able to determine how the run was configured, repeat its random-number sequence when possible, and judge whether differences across runs are consequential.
There are two useful goals. Computational repeatability means that the same data, model, software, settings, and random-number sequence reproduce the same calculations, often to numerical precision. Stochastic stability means that independently seeded runs give sufficiently similar estimates and conclusions for the intended use. A fixed seed can help with the first goal; it does not establish the second. Nor does agreement across a few runs prove that the algorithm found a global maximum.
For a basic description of the stochastic update, let \(S_k\) be the approximation to the sufficient statistics after iteration \(k\), and let \(\bar{s}_k\) be the statistic calculated from the latent-effect simulation at that iteration. A common update is
where \(\gamma_k\) is the gain. The M-step then updates the model parameters using the current statistic approximation. The random draws affect \(\bar{s}_k\); the gain schedule controls how strongly each new simulated result changes \(S_k\). Thus, both the random-number settings and the gain schedule are part of the algorithm being reproduced, not incidental details.
What to Record
A useful run record has enough detail to distinguish a true rerun from a run that merely uses the same nominal seed. Record the data version or a data fingerprint, model specification, parameterization, starting values, constraints, and any data transformations. Record software and relevant version information, including the implementation used for latent-effect simulation and numerical optimization. Differences in software versions, numerical libraries, or floating-point behavior can prevent bit-for-bit identity even when the statistical settings match.
For stochastic simulation, record the random-number generator (RNG) algorithm and the seed or saved RNG state. The seed initializes a sequence for a specified generator; a seed without the generator and its implementation may not identify the same sequence in another environment. Also record the simulation method, the number of latent-effect draws per iteration, any burn-in or transition settings, and whether draws are retained or discarded in particular phases. The number of draws matters because it influences the Monte Carlo variability of the quantities entering the stochastic update.
The gain schedule should be explicit: give the gain for each phase, the iteration at which a transition occurs, and any constants or exponent used in a formula. A common diminishing schedule has \(\gamma_k\) decrease after a finite initial phase. For the diminishing phase, a standard stochastic-approximation condition is \(0<\gamma_k\leq 1\), with
These conditions express that updates continue to have cumulative influence while their squared size becomes summable. They are not, by themselves, a guarantee that a particular finite run has converged or reached the desired solution. Documenting the exact schedule allows a reviewer to assess what was actually used rather than infer it from a label such as “standard SAEM settings.”
Finally, state the stopping rule in operational terms. Include the monitored quantities, tolerances, check frequency, number of consecutive checks required, minimum or maximum iteration count, and what happens if a maximum is reached without satisfying the criterion. If a run stops because an estimate becomes invalid or an optimizer fails, record that as a termination event, not as successful convergence.
Worked Examples
Worked Example: Same Seed, Different Run Records
Suppose two hypothetical analysts fit the same NLME model to the same cleaned data, each using seed 271828. Analyst A records the RNG algorithm and software version, uses 20 latent-effect draws per iteration, and runs 2,000 iterations. Analyst B records only the seed, uses 50 draws per iteration, and stops when a parameter-change criterion is met.
The seed alone does not make these the same computational run. The number of draws changes the simulation settings, and B's stopping rule may produce a different number of iterations. Even if the seed is identical, the random sequence can be consumed differently when the algorithms or software implementations differ. The two records therefore cannot support a claim of exact repeatability.
Now suppose Analyst A repeats the run with the same data, model, starting values, RNG algorithm, software environment, draw count, gain schedule, and fixed iteration count. If the implementation is deterministic given the RNG sequence, the repeat should reproduce the same trace, apart from possible small numerical differences. This tests repeatability of that configured run; it does not test whether its estimate is insensitive to the initial value or to Monte Carlo variation.
Worked Example: Reporting a Gain Update
In a hypothetical SAEM run, a single component of the sufficient-statistic vector has current value \(S_{k-1}=10.0\). The statistic from the simulated batch at iteration \(k\) is \(\bar{s}_k=14.0\), and the recorded gain is \(\gamma_k=0.25\). The updated value is
At the next iteration, suppose \(\bar{s}_{k+1}=9.0\) and \(\gamma_{k+1}=0.20\). Then
These calculations show why a run record needs both the simulated statistic and the gain sequence to explain the stochastic updates. Reporting only the final value \(10.6\) does not reveal whether the update was computed correctly or whether a different gain schedule was used. In an actual fit, the statistic is generally a vector, but the same update applies componentwise.
Worked Example: Applying a Fully Specified Stopping Rule
Consider a hypothetical rule that checks every 100 iterations after iteration 1,000. At each check, calculate the largest scaled change across monitored parameters and sufficient statistics over the last 100 iterations. Stop only if both quantities are below 0.01 at five consecutive checks. Let the scaled change be a maximum of absolute changes divided by prespecified scale constants, so parameters with different units can be assessed on a common relative scale.
At five successive checks, suppose the parameter changes are 0.008, 0.006, 0.009, 0.007, and 0.005, while the corresponding sufficient-statistic changes are 0.007, 0.009, 0.006, 0.008, and 0.004. Each value is strictly below 0.01, so both criteria pass at all five checks. If those checks occur at iterations 1,100, 1,200, 1,300, 1,400, and 1,500, the rule stops at iteration 1,500—not at the first passing check.
If the sufficient-statistic change at iteration 1,300 had instead been 0.012, the consecutive-check requirement would not be met there. The count of consecutive successful checks would reset, and five new consecutive passing checks would be needed. This illustrates why a report must state the check schedule and consecutive-check rule, not merely say “stopped when stable.”
Interpreting Differences Across Runs
A disciplined comparison separates three questions: Was the computational setup the same? Did the runs satisfy their stopping criteria? Are the resulting estimates meaningfully different on the scale of the model and its intended use? Matching seeds is useful for debugging an exact rerun, while different seeds are useful for examining sensitivity to the simulation path. These are different tasks and should be labeled as such.
For instance, suppose two hypothetical runs from different seeds estimate a clearance typical value of 4.82 L/h and 4.79 L/h. The absolute difference is \(4.82-4.79=0.03\) L/h, and relative to 4.82 L/h it is \(0.03/4.82\approx0.00622\), or about 0.622%. This small difference may be practically negligible, but its meaning depends on the model, precision needed, and uncertainty in the estimate. One pair of close estimates is not enough to establish general stability.
A useful comparison table can distinguish the role of each run:
| Run comparison | Main question | Interpretation |
|---|---|---|
| Same seed and unchanged setup | Can the configured computation be repeated? | Differences point to changed execution, numerical behavior, or an implementation issue. |
| Different seeds, same other settings | How sensitive are results to the stochastic path? | Spread reflects a combination of Monte Carlo variability and any remaining convergence or solution sensitivity. |
| Different initial values | Are results sensitive to initialization? | Distinct results may indicate different attraction regions or inadequate convergence; compare with the recorded stopping evidence. |
If independent-seed runs differ materially, inspect the traces and sufficient-statistic stabilization, not just the final estimates. Look for continued drift, noisy oscillation, boundary behavior, optimizer failures, or stopping at the iteration limit. A run that stopped under a threshold can still be inadequate if the threshold was poorly chosen or a monitored quantity was omitted. Conversely, modest differences between runs are expected in a stochastic algorithm and are not automatically evidence of a defect.
Common Pitfalls and Limits
- “We used the same seed, so the result must be reproducible.” A seed does not specify the RNG algorithm, software version, simulation settings, or the order in which random numbers are consumed. Record these details and the broader computational setup.
- “The estimates match, so convergence is proven.” Similar final estimates do not prove convergence, identify a global maximum, or show that all relevant diagnostics pass. Examine traces, stopping evidence, and sensitivity to initialization as appropriate.
- “The run converged” without defining the rule. A tolerance without the monitored quantity, scale, check frequency, and consecutive-check requirement cannot be independently applied. Report the full rule and the iteration at termination.
- Reporting only the final parameter estimates. Retain sufficient run information—such as parameter and statistic traces, termination status, and warnings—to diagnose how the output arose.
Reproducibility has limits. Different platforms or software versions may yield small numerical differences even with nominally identical settings, and a seed cannot make a poorly specified model or an inadequately converged run reliable. Documentation makes a result auditable; it does not replace statistical diagnostics or judgment about whether the estimates are useful.
Identify the data version, model, parameterization, starting values, and constraints.
Record the RNG algorithm and seed, simulation settings, and every phase of the gain schedule.
State monitored quantities, thresholds, check frequency, consecutive checks, and maximum-iteration behavior.
Save traces and termination details; use same-seed repeats for computational checks and different seeds to assess stochastic stability.
Check Your Understanding
Answer each question, then compare with the solutions below.
- Why does a seed by itself not fully specify a reproducible SAEM run?
- A gain update uses \(S_{k-1}=8\), \(\bar{s}_k=12\), and \(\gamma_k=0.30\). Calculate \(S_k\).
- A stopping rule requires four consecutive checks below 0.02. The observed values are 0.01, 0.015, 0.021, 0.008, 0.009, 0.012, and 0.019. Has the rule been met by the seventh check?
- What is the difference between using a repeated seed and using several different seeds?
Solutions. (1) A seed initializes a sequence for a particular RNG; the generator, software or implementation, simulation settings, gain schedule, inputs, and stopping rule also affect the run. (2) \(S_k=8+0.30(12-8)=8+1.2=9.2\). (3) No. The value 0.021 breaks the first sequence; the next four values are 0.008, 0.009, 0.012, and 0.019, all below 0.02, so the rule is met at the seventh check. (4) Repeating a seed under unchanged settings checks whether the same computation can be reproduced; different seeds probe sensitivity to stochastic simulation. Neither alone proves convergence.