Generating Replicate Data for Predictive Checks

From a Fitted Model to Replicate Outcomes

A predictive check asks whether data generated by a fitted model resemble the observed data in features relevant to the model’s intended use. Generating replicate outcomes is the step that turns the fitted model into hypothetical data: the model supplies the population parameters, the study supplies the design, and simulation supplies new subject-level and observation-level variation.

For a population predictive check, retain the study design that was actually observed: the subjects’ covariates, doses and dosing times, and observation times. Simulate new random effects and residual errors rather than copying the observed outcomes or treating estimated individual effects as known. This produces replicate outcomes at the same design points as the observed outcomes.

Definition: Let $i$ index subjects and $j$ index observations within subject $i$. A replicate outcome is $$ Y_{ij}^{\mathrm{rep}}=f(t_{ij};D_i,x_i,\hat{\theta},\eta_i^{\mathrm{rep}})+\epsilon_{ij}^{\mathrm{rep}}. $$ Here, $f$ is the structural model; $t_{ij}$ is the observation time; $D_i$ is subject $i$’s dosing history; $x_i$ denotes that subject’s covariates; $\hat{\theta}$ is the fitted population-parameter estimate; $\eta_i^{\mathrm{rep}}$ is a simulated random-effect vector; and $\epsilon_{ij}^{\mathrm{rep}}$ is simulated residual error. The precise way residual error enters the equation depends on the fitted error model.

In a standard population simulation, random effects are drawn from the fitted population distribution, commonly $\eta_i^{\mathrm{rep}}\sim N(0,\hat{\Omega})$, where $\hat{\Omega}$ is the estimated random-effects covariance matrix. Draw one vector for each subject, not one new vector for each observation. Residual errors are then generated according to the fitted residual model. For example, an additive model uses $Y_{ij}^{\mathrm{rep}}=m_{ij}^{\mathrm{rep}}+\epsilon_{ij}^{\mathrm{rep}}$, where $m_{ij}^{\mathrm{rep}}$ is the structural prediction and $\epsilon_{ij}^{\mathrm{rep}}\sim N(0,\hat{\sigma}^{2})$. A proportional model may use $Y_{ij}^{\mathrm{rep}}=m_{ij}^{\mathrm{rep}}(1+\epsilon_{ij}^{\mathrm{rep}})$, with $\epsilon_{ij}^{\mathrm{rep}}\sim N(0,\hat{\sigma}_{p}^{2})$.

These expressions assume the fitted model’s stated distributions and dependence structure are appropriate for the simulation. For instance, independent residual errors are appropriate only if that is the assumed residual structure; a model with serially correlated residuals requires errors simulated with that correlation. Likewise, the same random-effect draw applies to all records from a subject because it represents that subject’s latent deviation in the model.

What Is Held Fixed and What Is Simulated?

The distinction between design inputs and stochastic model components is central. For each replicate data set, hold the subjects, their covariates, doses, dose times, and observation times fixed when the target is prediction under the observed study design. Recalculate the structural predictions using those inputs and the simulated random effects. Then simulate residual errors and form the outcomes.

A common plug-in predictive check also holds $\hat{\theta}$ fixed. It therefore represents variability in outcomes conditional on the fitted population-parameter estimates, but it does not include uncertainty in those estimates. That may be the intended target. If the question requires parameter uncertainty as well, the simulation must draw population parameters from an uncertainty distribution or refit models to resampled data; simply drawing new random effects does not add parameter uncertainty.

Specify the simulation target.
Decide whether the replicate data should represent new subjects under the observed design, predictions conditional on individual information, or a target that also includes parameter uncertainty.

Keep the design inputs.
For the observed-design target, retain each subject’s covariates, dose amounts and times, and observation times.

Draw subject effects and compute predictions.
Draw one $\eta_i^{\mathrm{rep}}$ per subject from the fitted random-effects distribution and calculate the structural predictions at that subject’s scheduled times.

Draw residuals and form outcomes.
Simulate residual errors according to the fitted error model and combine them with structural predictions on the correct scale.

Check the simulated records.
Verify that records correspond to the intended design and that the generated values follow the model’s specified distributions and units.

Worked Simulations

Worked Example: One Subject With a Random Volume Effect

Hypothetical data and model. Suppose a subject receives a 100 mg intravenous bolus at time zero. Use a one-compartment model with clearance $CL=4$ L/h and typical volume $V_{\mathrm{typ}}=20$ L. Let $\log(V_i/V_{\mathrm{typ}})=\eta_{Vi}$, and suppose a simulated random effect is $\eta_{Vi}^{\mathrm{rep}}=\log(22/20)=0.09531$. Thus the simulated volume is $V_i^{\mathrm{rep}}=20e^{0.09531}=22$ L. The subject has observations scheduled at 1 h and 4 h. For illustration, use additive residual errors in mg/L, with simulated draws $-0.020$ and $+0.040$ mg/L.

The structural prediction after a bolus dose is $m(t)=(D/V_i^{\mathrm{rep}})\exp[-(CL/V_i^{\mathrm{rep}})t]$. The elimination rate is $CL/V_i^{\mathrm{rep}}=4/22=0.1818$ h$^{-1}$.

At 1 h, $$ m(1)=\frac{100}{22}\exp\left[-\frac{4}{22}(1)\right] =3.7898\ \text{mg/L}. $$ Adding the simulated residual gives $Y_{i1}^{\mathrm{rep}}=3.7898-0.020=3.7698$ mg/L, or 3.770 mg/L to three decimal places.

At 4 h, $$ m(4)=\frac{100}{22}\exp\left[-\frac{4}{22}(4)\right] =2.196478\ \text{mg/L}. $$ Thus $Y_{i2}^{\mathrm{rep}}=2.196478+0.040=2.236478$ mg/L, or 2.236 mg/L to three decimal places. The same simulated volume is used at both times; only the residual draw changes by observation.

Worked Example: Multiple Doses and Proportional Error

Hypothetical data and model. A subject receives 100 mg at time 0 and 50 mg at time 8 h. At 10 h, concentration is observed. Suppose this subject’s simulated volume is 25 L, clearance is 5 L/h, and the error model is proportional: $Y^{\mathrm{rep}}=m^{\mathrm{rep}}(1+\epsilon^{\mathrm{rep}})$, with a simulated residual draw $\epsilon^{\mathrm{rep}}=0.05$. Assume instantaneous bolus input and linear one-compartment elimination, so each dose contributes to the concentration remaining at the observation time.

The elimination rate is $k=CL/V=5/25=0.2$ h$^{-1}$. At 10 h, the first dose has been in the body for 10 h and the second for 2 h. Therefore, $$ m^{\mathrm{rep}}(10)=\frac{100}{25}e^{-0.2(10)}+\frac{50}{25}e^{-0.2(2)} =4(0.135335)+2(0.670320) =1.88198\ \text{mg/L}. $$

The proportional residual draw multiplies the prediction by $1.05$: $$ Y^{\mathrm{rep}}=1.88198(1+0.05)=1.97608\ \text{mg/L}. $$ The simulated outcome is approximately 1.976 mg/L. Applying the residual as an additive 0.05 mg/L would be incorrect: under this specified model, 0.05 is a dimensionless proportional error, not an amount in concentration units.

Worked Example: Replicates Are New Draws, Not Copies

Hypothetical setup. Consider one subject given a 100 mg bolus, with $CL=4$ L/h and $V_{\mathrm{typ}}=20$ L. The observations are scheduled at 1 h and 4 h. In one replicate, draw $\eta_V^{\mathrm{rep}}=0$, so $V=20$ L, and additive residuals $+0.100$ and $-0.100$ mg/L. The structural predictions are $$ m(1)=\frac{100}{20}e^{-(4/20)(1)}=4.09365\ \text{mg/L},\qquad m(4)=\frac{100}{20}e^{-(4/20)(4)}=2.24664\ \text{mg/L}. $$ The outcomes are $4.09365+0.100=4.19365$ mg/L and $2.24664-0.100=2.14664$ mg/L.

For a second replicate, draw a new random effect $\eta_V^{\mathrm{rep}}=\log(24/20)=0.18232$, giving $V=24$ L. With new residual draws $+0.050$ and $-0.050$ mg/L, the predictions become $$ m(1)=\frac{100}{24}e^{-(4/24)(1)}=3.5270\ \text{mg/L},\qquad m(4)=\frac{100}{24}e^{-(4/24)(4)}=2.1392\ \text{mg/L}. $$ The replicate outcomes are approximately 3.577 mg/L and 2.089 mg/L. Both records within each replicate share that replicate’s volume, while the second replicate receives a new volume draw. The scheduled times remain 1 h and 4 h in both.

Why the Simulation Represents the Predictive Distribution

With fitted population parameters treated as fixed, the model’s predictive distribution for a subject’s outcome averages over possible random effects and residual errors. In notation, for design information $D_i,x_i,t_{ij}$, this is the distribution obtained by integrating the conditional outcome distribution over the random-effects distribution: $$ p(y_{ij}\mid D_i,x_i,t_{ij},\hat{\theta},\hat{\Omega}) =\int p(y_{ij}\mid \eta_i,D_i,x_i,t_{ij},\hat{\theta})\,p(\eta_i\mid\hat{\Omega})\,d\eta_i. $$

This integral is often difficult to evaluate directly for nonlinear models. Simulation approximates it by repeatedly drawing $\eta_i^{\mathrm{rep}}$, then drawing residual error conditional on the resulting prediction. The procedure is justified because the draws follow the model’s assumed generative sequence: first subject effects, then model predictions given those effects and the design, then observation errors given the predictions.

Checks, Alternatives, and Limitations

Before using replicate outcomes in a predictive check, confirm that the simulation matches the intended target. Useful checks include verifying the number and identity of records per subject, confirming doses and times were applied in the correct sequence, and checking that one random-effect draw is shared across all of a subject’s observations. Inspect simulated values for impossible values or unit inconsistencies. For example, an additive normal error can generate negative concentrations; if the model is defined on the original concentration scale, that possibility should be assessed rather than silently truncating values, which changes the model’s implied distribution.

A population simulation that draws new random effects is not the same as a simulation conditional on individual empirical Bayes estimates. Conditional simulations hold subject-specific information fixed or use it to inform subject effects, and may be appropriate for questions about individual predictions. They are not interchangeable with population predictive simulations: shrinkage of individual estimates can make the simulated between-subject variation too small if those estimates are used as though they were the population distribution.

The main alternative to a plug-in simulation is to include uncertainty in population parameter estimates, for example by drawing parameter sets from an estimated sampling or posterior distribution and simulating outcomes for each set. This addresses a broader uncertainty target, but requires a defensible parameter-uncertainty distribution. Neither approach repairs a misspecified structural model, random-effects distribution, residual model, or dosing record. Replicate data reflect the assumptions supplied to the simulator, not the assumptions one intended to supply.

Common pitfall: Simulating one new random effect for every observation treats repeated measurements from the same subject as if they belonged to different subjects. The correction is to draw $\eta_i^{\mathrm{rep}}$ once per subject and reuse it across that subject’s records; draw observation-level residual errors separately according to the fitted residual structure.

A further pitfall is to simulate at convenient, evenly spaced times rather than the actual observation times. That changes the design and can change the distribution of summaries, particularly when sampling is sparse or uneven. For a check aimed at the observed study, retain its actual times and dosing histories. Alternative designs can be simulated when they are themselves the explicit target.

Check Your Understanding

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

In a standard population replicate, which quantities are ordinarily held fixed and which are simulated?
A subject has five observations and a random intercept. How many random-intercept draws should be used for that subject in one replicate, and why?
For a proportional error model, the prediction is 8 mg/L and the simulated dimensionless residual is $-0.10$. What replicate outcome results?
Why does using individual empirical Bayes estimates in place of newly drawn random effects change the predictive target?
Do plug-in simulations with fixed fitted population parameters include uncertainty in those parameter estimates?

Solutions. (1) For an observed-design check, hold the subjects, covariates, doses and dose times, and observation times fixed; simulate random effects and residual errors, while usually holding fitted population parameters fixed. (2) Draw one random intercept and reuse it for all five observations, because it represents a subject-level deviation. (3) $8(1-0.10)=7.2$ mg/L. (4) Those estimates are conditioned on each subject’s data and can be shrunk; they do not represent new draws from the fitted population distribution. (5) No. Fixed-parameter simulations include outcome variation under the fitted model, not uncertainty in the estimated population parameters.