Why Starting Values Matter in SAEM
SAEM estimates population parameters when individual random effects are latent. The observed-data likelihood therefore involves integrating over effects that are not directly observed. In nonlinear models, that likelihood can have a complicated surface: distinct parameter combinations may fit the data reasonably well, and some may correspond to different local maxima. An algorithm can settle near one such solution without finding the highest likelihood available.
Starting values matter in two related but distinct ways. Initial population parameters set the algorithm's current location in parameter space. Initial latent-effect states set the starting point of the simulation used to approximate the latent quantities. Either can affect the path taken during estimation. A stable result from one run is not, by itself, evidence that the same solution would be found from other plausible starts.
At a smooth interior local maximum, the likelihood gradient is zero and the Hessian is negative semidefinite; for a strict, well-identified maximum, it is typically negative definite. Constraints change the boundary conditions: a maximum at a boundary need not have zero gradient in every unconstrained direction. These characterizations describe candidate maxima, not a guarantee that a numerical algorithm has found one.
How Parameters and Latent States Enter
Write the observed data for subject \(i\) as \(y_i\), the subject's latent random effects as \(\eta_i\), and the population parameters as \(\boldsymbol{\theta}\). The observed-data likelihood has the form
SAEM alternates simulation of latent effects with stochastic updating of sufficient statistics and a parameter update. The parameter start \(\boldsymbol{\theta}^{(0)}\) affects the conditional distribution used to simulate the effects; the initial states \(\eta_i^{(0)}\) affect where those simulations begin. They are not the same choice. A sound latent-state initialization does not repair poor population-parameter values, and sensible population parameters do not guarantee that latent-effect simulation explores the relevant conditional distribution.
If latent-effect simulation mixes adequately and the algorithm's regularity and stochastic-approximation conditions hold, the influence of its initial states can diminish. But slow mixing, weakly identified parameters, boundaries, or a poorly explored likelihood surface can leave practical dependence on the start. SAEM's stochastic updates also mean that the likelihood need not increase at every individual iteration. A settled trace from one run is therefore a convergence check, not a global-optimality certificate.
Worked Example: Translating a Population Start Into Individual Clearances
Consider a hypothetical clearance model with \(CL_i=TVCL\exp(\eta_i)\), where \(TVCL\) is the typical clearance in L/h and \(\eta_i\) is a dimensionless subject effect. Suppose a run starts with \(TVCL=6.0\) L/h. Compare the subject-specific initial clearance for \(\eta_i=0\) and \(\eta_i=0.30\).
For \(\eta_i=0\), \(CL_i=6.0\exp(0)=6.0\) L/h. For \(\eta_i=0.30\),
The population value is the same in both initializations, but the latent-effect state implies a different starting individual clearance. This can change the initial conditional distributions and thus the early SAEM trajectory. The calculation does not establish that either state is a good estimate: it only makes explicit what the chosen start implies under this model.
Why Different Runs Can Find Different Solutions
The EM idea helps explain the mechanism. At a current parameter value, an EM update maximizes an expectation of the complete-data log-likelihood, conditional on the observations. In SAEM, simulation and stochastic averaging approximate the needed latent-effect expectations. If the likelihood surface has separated attractive regions, updates initialized in different regions may approach different stationary solutions. With SAEM, finite simulation and gain schedules add stochastic variation to this path.
A useful local justification comes from a smooth log-likelihood near a candidate \(\boldsymbol{\theta}^{*}\). For a small displacement \(\boldsymbol{h}\), a second-order expansion gives
At an interior stationary point, \(\nabla\ell(\boldsymbol{\theta}^{*})=0\). If the Hessian \(H\) is negative definite, then \(\boldsymbol{h}^{T}H\boldsymbol{h}<0\) for every nonzero sufficiently small displacement, so the likelihood decreases locally. This establishes a strict local maximum, not a global one: it says nothing about a distant parameter value with a higher likelihood.
For SAEM, the practical response is not to infer global optimality from convergence alone. Use plausible, meaningfully varied initial values, examine whether solutions cluster, and compare final candidates using the same observed-data likelihood evaluation and the same data and model. Parameter plausibility and boundary behavior matter alongside the likelihood.
Worked Example: Comparing Three Converged Starts
Suppose three hypothetical SAEM runs use the same nonlinear clearance model and data. Their converged estimates and log-likelihoods, evaluated with the same high-accuracy likelihood method, are shown below. All runs satisfy the model's parameter constraints and passed the same convergence checks.
| Run | Initial \(TVCL\) (L/h) | Final \(TVCL\) (L/h) | Final \(\omega_{CL}\) | Final \(\ell\) |
|---|---|---|---|---|
| A | 5.0 | 6.2 | 0.22 | -128.4 |
| B | 10.0 | 8.0 | 0.35 | -130.1 |
| C | 6.5 | 6.3 | 0.23 | -128.5 |
Runs A and C yield similar estimates and log-likelihoods: their log-likelihood difference is \((-128.4)-(-128.5)=0.1\). Run B gives a different parameter combination and a lower log-likelihood. The likelihood ratio of A to B is
Thus the fitted likelihood at A is about 5.474 times the fitted likelihood at B, under this common evaluation. This is a comparison of candidate fits, not a likelihood-ratio hypothesis test: the runs fit the same model, and the statistic does not have a stated reference distribution for testing these two local solutions. A and C provide evidence of a repeatable solution region; B signals that the result may depend on initialization. More starts and checks are needed before treating A as the preferred solution, and these results alone do not prove it is the global maximum.
Choosing and Assessing Multiple Starts
A useful multiple-start exercise varies both population parameters and, when relevant, initial latent states. Choose starts that are scientifically plausible and cover distinct regions of the admissible space; arbitrary extreme values can instead create invalid or numerically uninformative runs. Keep the data, model, constraints, and likelihood-comparison method fixed. Record whether each run converged, its final parameter estimates, its objective value, and any boundary or numerical warnings.
Choose several plausible population-parameter vectors, including meaningfully different values for weakly identified or nonlinear parameters. If latent states are initialized explicitly, document those choices separately.
Apply the same estimation settings and stopping criteria. Inspect parameter traces and sufficient-statistic stabilization, and flag runs with nonconvergence, invalid covariance matrices, or boundary estimates.
Evaluate final candidates with a consistent observed-data likelihood method. Compare estimates as well as likelihoods; close objective values can accompany materially different parameters when the model is weakly identified.
State how many starts were used, how many reached acceptable solutions, whether solutions clustered, and why a candidate was selected. Do not claim proof of a global maximum from a finite set of starts.
Worked Example: Separating a Latent-State Start From a Parameter Start
Consider a hypothetical random-intercept model \(y_{ij}=\theta+\eta_i+\epsilon_{ij}\), with \(\eta_i\sim N(0,\omega^2)\) and independent \(\epsilon_{ij}\sim N(0,\sigma^2)\). For one subject, suppose the three observations are 11, 9, and 12, and use starting values \(\theta=10\), \(\omega^2=4\), and \(\sigma^2=1\). Find the conditional mode of \(\eta_i\) as a possible initial latent state.
The conditional density is proportional to the product of the observation density and the random-effect density. Ignoring constants, its log is
Differentiate with respect to \(\eta_i\), set the derivative to zero, and solve:
Here the residual sum before the random intercept is \((11-10)+(9-10)+(12-10)=2\). Therefore
Because the conditional log density is a quadratic with negative second derivative \(-3/\sigma^2-1/\omega^2=-3.25\), this stationary point is its unique conditional maximum for these fixed parameter values. Initializing \(\eta_i\) at 0 or at 0.6154 changes the latent-state start, not \(\theta\), \(\omega^2\), or \(\sigma^2\). The conditional mode is a possible convenient initialization, not a substitute for simulating the conditional distribution during SAEM.
Diagnostics, Alternatives, and Limitations
Signs that starting values may be influencing the result include materially different final parameter vectors across runs, distinct likelihood levels, repeated convergence to a boundary, or traces that appear stable only for some starts. Also investigate candidates with near-identical likelihoods but markedly different estimates: this can indicate a flat or weakly identified direction rather than two sharply separated maxima. A run that stops without stabilized parameter and sufficient-statistic traces should not be counted as a successful converged solution merely because it returned estimates.
Multiple starts are a practical search strategy, not a proof procedure. They can reveal competing solutions found by the chosen algorithm, but may miss other regions. Another approach is to use deterministic optimization or a different likelihood approximation as a check, where suitable. Such methods have their own numerical and modeling limitations; agreement is reassuring but still not a mathematical guarantee of global optimality. SAEM is useful for the latent-variable estimation problem, while the multiple-start comparison addresses sensitivity to initialization. These are complementary, not competing, roles.
A common mistake is to treat the run with the largest reported objective as automatically correct, even if its likelihood was evaluated differently, its parameters violate constraints, or the run failed convergence checks. The correction is to compare only admissible, adequately converged candidates under a consistent likelihood evaluation, inspect their parameter estimates and diagnostics, and describe the conclusion as evidence from the starts examined—not proof that no better solution exists.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- In a model with \(CL_i=TVCL\exp(\eta_i)\), what individual clearance is implied by \(TVCL=5.0\) L/h and \(\eta_i=0.20\)?
- Two valid runs have log-likelihoods \(-80.0\) and \(-81.0\). What is the likelihood ratio of the first fit to the second, and what does it mean?
- Why does a stable parameter trace from a single SAEM run not establish that the global maximum was found?
- Give one possible interpretation if multiple runs have nearly equal likelihoods but noticeably different parameter estimates.
Solutions. 1. \(CL_i=5.0\exp(0.20)=5.0(1.2214028\ldots)\approx 6.107\) L/h. 2. The ratio is \(\exp\{(-80.0)-(-81.0)\}=\exp(1)\approx 2.718\); under the same likelihood evaluation, the first fitted likelihood is about 2.718 times the second. This is not by itself a hypothesis test. 3. A stable trace indicates settling in that run, but other regions of the likelihood surface may contain a higher maximum. 4. The likelihood may be relatively flat in some parameter directions, suggesting weak identification; inspect uncertainty and model diagnostics rather than concluding that the estimates are interchangeable.