A Stepwise Algorithm for Simulating From a Fitted NLME Model

Why Simulation Needs an Order

An NLME model describes a linked set of quantities, not a collection of unrelated random draws. Population parameters determine how covariates and subject effects contribute to predictions; the prediction, in turn, determines the distribution of an observation when residual error depends on its magnitude. A simulation is coherent only when each quantity is generated from the distribution specified for it by the model.

The appropriate sequence is therefore more than a programming convenience. It encodes the assumptions behind the simulated data. For example, using one population-parameter draw for an entire simulated dataset represents uncertainty about the population model, whereas independently redrawing those parameters for each subject represents a different, generally unintended, model.

Definition: A simulation path is one set of generated values of the quantities required to produce a dataset: population parameters \(\theta\), covariates \(x_i\), subject random effects \(\eta_i\), model predictions, and residual errors \(\epsilon_{ij}\). A simulated observation \(Y_{ij}\) is generated from its conditional distribution given the relevant parameters, covariates, random effects, and design.

The Generative Model and Its Factorization

Use \(i\) to index subjects and \(j\) to index records. Let \(\theta\) denote the population parameters, \(x_i\) the covariates for subject \(i\), \(\eta_i\) that subject’s random effects, and \(t_{ij}\) the design information for record \(j\), such as its sampling time. Write the model prediction as \(f(t_{ij},x_i,\theta,\eta_i)\). A common observation model is

Y_{ij}=f(t_{ij},x_i,\theta,\eta_i)+\epsilon_{ij}.

This equation separates the model prediction from the residual error. The prediction is the model’s conditional mean when the residual has mean zero; it is not yet the simulated observation. The distribution of \(\epsilon_{ij}\) may depend on the prediction or other model quantities, so it must be specified conditionally rather than assumed to be the same in every setting.

For a population simulation that draws uncertain parameters, one useful general factorization is

p(\theta,x,\eta,Y) = p(\theta)\, p(x\mid\theta)\, p(\eta\mid\theta,x)\, p(Y\mid\theta,x,\eta).

Here \(x\), \(\eta\), and \(Y\) stand for the covariates, random effects, and observations over the simulated dataset. This factorization follows the chain rule for a joint distribution; it does not assert that covariates are independent of parameters or that random effects are independent of covariates. If a simulation conditions on known covariates, \(p(x\mid\theta)\) is a point mass at those values. If covariates are independent of \(\theta\), then \(p(x\mid\theta)=p(x)\), but that is an additional assumption, not an automatic feature of NLME models.

Some models make further conditional-independence assumptions. For instance, the usual subject-level random-effects model may specify that \(\eta_i\) are independent across subjects conditional on the population parameters and covariates. A residual model may specify conditional independence across records. Use those assumptions only when they are part of the fitted model and simulation target. When a model instead gives a joint distribution for dependent quantities, draw them jointly or use the relevant conditional distributions.

A Stepwise Simulation Algorithm

Before drawing anything, state the target of the simulation. Decide whether \(\theta\) is fixed at its fitted estimate or drawn from a distribution representing parameter uncertainty; decide whether covariates are fixed or sampled; and specify which subjects and records the dataset contains. The algorithm below expresses the generative order while allowing a step to be omitted when its value is conditioned on.

Choose the population-parameter value.
Either set \(\theta\) to the fitted estimates, or draw one \(\theta\) from the chosen uncertainty distribution for the simulated dataset. Keep the interpretation of that choice explicit.

Set or draw covariates and design.
Use known covariates and sampling times when conditioning on them. Otherwise draw covariates from the specified \(p(x\mid\theta)\), including any dependence on \(\theta\).

Draw subject random effects.
For each subject, draw \(\eta_i\) from the model’s specified distribution conditional on the quantities it depends on, such as \(\theta\) and \(x_i\). Generate a subject’s random effect once for the subject-level records it governs.

Calculate model predictions.
Evaluate \(f(t_{ij},x_i,\theta,\eta_i)\) for each required record using the values already assigned to that dataset and subject.

Generate residual errors and observations.
Draw residual errors from the specified conditional observation model, then combine them with the predictions to produce \(Y_{ij}\).

The sequence is a practical way to respect the factorization, not a claim that all random variables must be generated one at a time. A joint draw is appropriate when the model specifies a joint distribution. The essential rule is that a quantity must be drawn from the correct distribution conditional on the quantities on which it depends.

Worked Paths Through the Algorithm

Worked Example: From a Population Model to One Observation

Suppose a hypothetical one-compartment model after an intravenous bolus has prediction

f(t)=\frac{D}{V_i}\exp\left(-\frac{CL_i}{V_i}t\right), \qquad CL_i=CL\exp(\eta_{CL,i}),\quad V_i=V\exp(\eta_{V,i}).

For this example, condition on a dose \(D=400\ \mathrm{mg}\) and a time \(t=4\ \mathrm{h}\). Set the fitted population values to \(CL=3.2727\ \mathrm{L/h}\) and \(V=40.0\ \mathrm{L}\). Consider one hypothetical random-effects draw, \(\eta_{CL,i}=\eta_{V,i}=\ln(1.1)\), from the model’s assumed joint random-effects distribution. The draw is illustrative; the simulation must use the actual distribution specified for the application.

First calculate the individual parameters:

CL_i=3.2727\ \mathrm{L/h}\times 1.1 \approx 3.6000\ \mathrm{L/h}, \qquad V_i=40.0\ \mathrm{L}\times 1.1 =44.0\ \mathrm{L}.

Next calculate the prediction, using hours and litres consistently:

\begin{aligned} \mu_i(4) &= \frac{400\ \mathrm{mg}}{44.0\ \mathrm{L}} \exp\left(-\frac{3.6000\ \mathrm{L/h}}{44.0\ \mathrm{L}}(4\ \mathrm{h})\right)\\ &= 9.0909\ \mathrm{mg/L}\times e^{-0.32727}\\ &\approx 6.554\ \mathrm{mg/L}. \end{aligned}

Assume additive residual error \(\epsilon_{ij}\sim N(0,(0.200\ \mathrm{mg/L})^2)\), conditional on this prediction. For one hypothetical residual draw of \(-0.200\ \mathrm{mg/L}\), the observation is

Y_{ij}=6.55352\ \mathrm{mg/L}-0.200\ \mathrm{mg/L} =6.35352\ \mathrm{mg/L} \approx 6.354\ \mathrm{mg/L}.

The unrounded prediction is retained in the addition, so the reported observation is consistent with the calculation. The resulting observation is one possible draw, not a replacement for the prediction or a claim about any real study.

Worked Example: Covariates That Depend on the Parameter Draw

Suppose a hypothetical population simulation has two possible parameter states, \(A\) and \(B\), each with probability \(0.5\). Let \(X=1\) indicate membership in a covariate group. Assume \(P(X=1\mid A)=0.8\) and \(P(X=1\mid B)=0.2\). These probabilities explicitly make the covariate distribution depend on the parameter state.

The joint probability of drawing state \(A\) and \(X=1\) is

P(A,X=1)=P(A)P(X=1\mid A)=0.5(0.8)=0.40.

The marginal probability of \(X=1\) is \(0.5(0.8)+0.5(0.2)=0.50\). If one incorrectly treated \(X\) as independent of the parameter state, the joint probability would be \(P(A)P(X=1)=0.5(0.5)=0.25\), not \(0.40\). Thus, drawing from the marginal covariate distribution after choosing a parameter state would fail to reproduce the stated joint model.

For a concrete path, suppose state \(A\) is drawn, then \(X=1\), then \(\eta=-0.4\ \mathrm{mg/L}\). Let the prediction be \(f=10+2X+\eta\ \mathrm{mg/L}\), with the coefficient \(2\) also in \(\mathrm{mg/L}\), and let a residual draw be \(+0.3\ \mathrm{mg/L}\). The prediction is \(10+2(1)-0.4=11.6\ \mathrm{mg/L}\), and the observation is \(11.6+0.3=11.9\ \mathrm{mg/L}\). The covariate draw follows the selected state, and the observation follows the resulting prediction.

Worked Example: One Population Draw Shared Across a Dataset

Consider a hypothetical two-subject dataset from the model \(f_i=\theta_0+\theta_1x_i+\eta_i\), where \(x_i\) is dimensionless and \(f_i,\eta_i,\theta_0,\theta_1\) are in \(\mathrm{mg/L}\). Suppose one parameter draw for the dataset gives \(\theta_0=8.0\ \mathrm{mg/L}\) and \(\theta_1=1.5\ \mathrm{mg/L}\). Condition on \(x_1=0\) and \(x_2=2\). For illustration, use subject-effect draws \(\eta_1=0.4\ \mathrm{mg/L}\) and \(\eta_2=-0.2\ \mathrm{mg/L}\), followed by residual draws \(\epsilon_1=0.1\ \mathrm{mg/L}\) and \(\epsilon_2=-0.3\ \mathrm{mg/L}\).

For subject 1, the prediction is \(8.0+1.5(0)+0.4=8.4\ \mathrm{mg/L}\), so \(Y_1=8.4+0.1=8.5\ \mathrm{mg/L}\). For subject 2, the prediction is \(8.0+1.5(2)-0.2=10.8\ \mathrm{mg/L}\), so \(Y_2=10.8-0.3=10.5\ \mathrm{mg/L}\). The same population-parameter draw is used for both subjects, while each subject has a separate random-effect draw. Redrawing \(\theta_0\) and \(\theta_1\) independently for each subject would describe a different population-level process.

Why the Order Works

The justification is the chain rule for the joint distribution. For any ordered set of random quantities, the joint density or probability can be written as a product of successive conditional distributions. Applying that rule in the order \(\theta,x,\eta,Y\) gives

p(\theta,x,\eta,Y) = p(\theta)\, p(x\mid\theta)\, p(\eta\mid\theta,x)\, p(Y\mid\theta,x,\eta).

A simulation that first draws \(\theta\), then draws \(x\) from \(p(x\mid\theta)\), then \(\eta\) from \(p(\eta\mid\theta,x)\), and finally \(Y\) from \(p(Y\mid\theta,x,\eta)\) therefore generates from this joint distribution. If a specified conditional distribution does not depend on one of its conditioning variables, that dependence can be omitted. For example, if \(x\) is independent of \(\theta\), then \(p(x\mid\theta)=p(x)\). The factorization itself does not establish that independence; it accommodates dependence when present.

For a fitted-model simulation with parameters fixed at their estimates, take \(p(\theta)\) to be a point mass at those estimates. For a simulation that includes parameter uncertainty, \(p(\theta)\) is the chosen distribution of parameter draws. The remaining steps still follow the conditional model. This distinction separates uncertainty about the population parameter values from subject-to-subject variability and residual variability.

Diagnostics, Alternatives, and Limitations

A useful audit records, for a few simulated paths, the population draw, covariates, random effects, predictions, residuals, and final observations. Check that the arithmetic reproduces the observations and that each value has the intended units. Also examine whether simulated values are plausible under the model: unexpected scales, impossible values, or patterns inconsistent with the stated conditional distributions can reveal an incorrect transformation, an omitted dependency, or a coding error.

Check dependencies. If covariates depend on \(\theta\), sampling them from an unconditional distribution can give the wrong joint population. Use the conditional distribution or explicitly condition on known covariates.
Check the level of each draw. Population parameters, subject effects, and residual errors do not automatically vary at the same level. Generating a new population-parameter draw for each subject changes the model being simulated.
Check the prediction before adding error. A residual draw changes an observation, not the model prediction. Do not feed a residual-contaminated observation into the prediction equation unless the model specifically requires that construction.
Check the residual model. Additive, proportional, and combined error models imply different conditional distributions. Drawing residuals from a constant-variance normal distribution is incorrect when the specified residual model says otherwise.

The main alternative is a plug-in simulation: set \(\theta\) to its fitted estimate and omit parameter uncertainty. This is appropriate when the target is simulation conditional on those estimates, but it does not represent uncertainty in them. A simulation using parameter draws instead incorporates the uncertainty distribution selected for the task; the quality of that result depends on the appropriateness of that distribution. Likewise, conditioning on observed covariates answers a different question from sampling covariates for a population. Neither choice is universally preferable: the simulation design must match the scientific target.

The algorithm also inherits the limitations of the fitted model and the chosen input distributions. It cannot correct model misspecification, an inadequate covariate distribution, or an unsuitable residual model. A syntactically correct sequence of random draws is not evidence that the model is scientifically adequate; it only ensures that the simulation follows the assumptions supplied to it.

Common pitfall: A student draws a parameter value, samples covariates from their marginal distribution, draws random effects, and adds residual error, then claims the procedure follows the model because each component has the right marginal distribution. That reasoning is wrong when the components are dependent: matching each marginal distribution does not guarantee the correct joint distribution. Draw each quantity from its specified conditional distribution, or draw dependent quantities jointly.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

Write the general factorization for a simulation that allows covariates to depend on population parameters.
In a one-subject model, the prediction is \(12.0\ \mathrm{mg/L}\) and a residual draw is \(-0.5\ \mathrm{mg/L}\). What is the simulated observation?
When is it valid to replace \(p(x\mid\theta)\) by \(p(x)\)?
For a two-subject dataset, which quantity is ordinarily shared when one population-parameter draw is used for the dataset: \(\theta\) or each subject’s \(\eta_i\)?
Why can correct marginal distributions for covariates and parameters still produce an incorrect simulation?

Solutions. (1) \(p(\theta,x,\eta,Y)=p(\theta)p(x\mid\theta)p(\eta\mid\theta,x)p(Y\mid\theta,x,\eta)\). (2) \(Y=12.0-0.5=11.5\ \mathrm{mg/L}\). (3) When covariates are independent of \(\theta\) under the simulation model, or when the relevant conditioning is otherwise immaterial. (4) The population draw \(\theta\) is shared; subject effects are drawn according to their subject-level model. (5) Because dependencies determine joint probabilities, and the wrong dependence structure changes the joint distribution even if the separate marginal distributions are correct.