Generating Replicate Data From a Population PK Model

What a Replicate Dataset Represents

A predictive check compares observed data with data the fitted model could plausibly generate. A replicate dataset is one simulated realization of the model under a specified design. It is not a copy of the observed concentrations, nor is it a set of predictions with residual error removed. It includes the sources of variation that belong to the predictive target: potentially population-parameter uncertainty, between-subject variability, residual error, and variation in doses or observation times.

The choice of what to simulate depends on the question. For a population-level check, a replicate typically represents a new set of individuals from the modeled population, often with the observed dosing and sampling design held fixed. For a conditional predictive check, the simulation instead conditions on subject-specific information, such as estimated individual effects, according to the definition of the check. These targets were distinguished in the earlier material on population-level and conditional predictive checks; the essential point here is to make the conditioning explicit before generating data.

Definition: A replicate dataset is a simulated dataset generated from the specified population PK model and design, using a declared set of fixed and simulated quantities. A valid replicate must preserve the dependence structure in the model: population parameters are shared at the appropriate level, individual random effects vary by individual, and residual errors vary at the observation level.

A convenient general representation is \(Y_{ij}^{\mathrm{rep}}=h(\theta,\eta_i,d_i,t_{ij})+\text{residual error}\), where \(i\) indexes individuals, \(j\) indexes observations, \(\theta\) denotes population parameters, \(\eta_i\) denotes individual random effects, \(d_i\) represents dose history, and \(t_{ij}\) is an observation time. The function \(h\) is the model-predicted concentration. The exact form of \(h\) and the error model must be those of the fitted model, not an alternative chosen for convenience during simulation.

Specify the Generative Model First

Consider a one-compartment model after an intravenous bolus. For a single dose \(D_i\) at time zero, the concentration predicted for individual \(i\) at time \(t\) is

C_i(t)=\frac{D_i}{V_i}\exp\left(-\frac{CL_i}{V_i}t\right).

Here \(CL_i\) is clearance in L/h, \(V_i\) is volume in L, \(D_i\) is dose in mg, and \(C_i(t)\) is in mg/L. One common lognormal random-effects specification is

CL_i=CL_{\mathrm{pop}}\exp(\eta_{CL,i}),\qquad V_i=V_{\mathrm{pop}}\exp(\eta_{V,i}),\qquad \begin{pmatrix}\eta_{CL,i}\\\eta_{V,i}\end{pmatrix}\sim N(\mathbf{0},\Omega).

The matrix \(\Omega\) contains the variances and covariance of the random effects. Its diagonal elements must be nonnegative and it must be positive semidefinite. A lognormal residual model can be written as

Y_{ij}=C_i(t_{ij})\exp(\epsilon_{ij}),\qquad \epsilon_{ij}\sim N(0,\sigma^2).

Under this specification, residual error is multiplicative and conditionally independent across observations given the model quantities. The simulated observation is positive. Other fitted models may use additive, combined, or different residual-error structures; use the fitted structure and its parameterization exactly. In particular, do not add error twice if the model software's prediction already includes a specified observation model.

Before simulating, specify whether \(\theta\) is fixed at its estimate or drawn from a parameter-uncertainty distribution. A plug-in simulation uses the fitted population estimates for every replicate. A simulation that propagates parameter uncertainty draws one plausible parameter vector per replicate, for example from a suitable posterior distribution or an approximate sampling distribution on a scale that respects parameter constraints. That vector is shared by all individuals in that replicate. Drawing a new population clearance independently for every individual would incorrectly turn population-parameter uncertainty into between-subject variability.

From the Model to a Replicate

The simulation is a sequence of conditional draws. Population parameters are selected first; then individual effects are generated (or the conditional effects are specified for a conditional check); concentrations are computed from the dose histories and observation times; and residual errors are added according to the observation model. If the target includes a stochastic design, doses and times are generated from its design mechanism as well. If the target conditions on the observed design, use the observed doses and times instead.

Set the predictive target.
State which population, individuals, covariates, doses, and sampling times the replicate represents. Identify what is fixed and what is simulated.

Draw population parameters if required.
Use one parameter vector for the replicate. If using a distribution for parameter uncertainty, draw on a scale and from a distribution that preserve parameter constraints and estimated dependence.

Generate individual effects.
For a population-level replicate, draw one random-effects vector per individual from the specified distribution using the selected \(\Omega\). Preserve modeled correlations.

Construct each concentration prediction.
Apply the PK model to each individual's parameters, covariates, dose history, and observation times.

Generate observations.
Draw residual errors at the observation level and apply the fitted residual model. Apply any modeled observation or sampling mechanism that belongs to the predictive target.

The justification is the law of total probability. If \(\theta\) is uncertain, \(\eta\) denotes the collection of individual effects, \(D,T\) denote design quantities, and \(Y^{\mathrm{rep}}\) denotes replicated observations, then the predictive distribution is obtained by averaging the conditional data distribution over the quantities drawn:

p(Y^{\mathrm{rep}}\mid \text{data}) = \int p(Y^{\mathrm{rep}}\mid\theta,\eta,D,T)\, p(\eta,\theta,D,T\mid\text{data})\, d\eta\,d\theta\,dD\,dT.

For a plug-in check, \(p(\theta\mid\text{data})\) is replaced by a point mass at the fitted estimate; for a fixed design, the design distribution is a point mass at the specified doses and times. Simulation approximates this integral by drawing the latent quantities and then drawing observations conditionally. For new individuals in a population-level check, effects are instead drawn from the population random-effects distribution \(p(\eta\mid\theta,D,T)\); for a conditional check, use the data-informed distribution \(p(\eta\mid\theta,D,T,\text{data})\). This distinction also explains why the order and level of the draws matter.

Worked Example: Individual Effects and Residual Error

Suppose, hypothetically, that a plug-in replicate uses \(CL_{\mathrm{pop}}=5.00\) L/h and \(V_{\mathrm{pop}}=25.0\) L. For one individual, the simulated random effects are \(\eta_{CL}=\log(1.10)\) and \(\eta_V=\log(0.80)\). The individual receives a 200 mg IV bolus, and a concentration is simulated at 2.00 h. Let the residual standard deviation on the log scale be \(\sigma=0.100\), and let the simulated standardized residual draw be \(z=-0.500\).

The individual parameters are \(CL_i=5.00\exp(\log(1.10))=5.50\) L/h and \(V_i=25.0\exp(\log(0.80))=20.0\) L. Therefore the elimination rate constant is \(CL_i/V_i=5.50/20.0=0.275\) h\(^{-1}\), and

C_i(2.00) =\frac{200}{20.0}\exp(-0.275\times2.00) =10.0\exp(-0.550) =5.76950\ \text{mg/L}.

The residual draw is \(\epsilon=\sigma z=0.100(-0.500)=-0.0500\). Applying the lognormal observation model gives

Y^{\mathrm{rep}}=5.76950\exp(-0.0500)=5.48812\ \text{mg/L}.

This is one simulated observation in one hypothetical replicate. The individual effect changes the PK parameters before the concentration is calculated; residual error is then applied to the predicted concentration. The calculation does not imply that every individual in the replicate has the same random effects or that every observation has the same residual draw.

Including Parameter Uncertainty

Parameter uncertainty is a separate source of variation from individual random effects. If it is part of the predictive target, use a joint distribution for the fitted population parameters. One option is a posterior distribution; another is an approximate sampling distribution from estimation, when that approximation is appropriate. Positive parameters are commonly represented on a log scale, and a covariance structure should be retained when multiple estimates are drawn together. A normal approximation on the original scale can generate impossible negative values, so it should not be used mechanically.

Worked Example: A Shared Population-Parameter Draw

For illustration, suppose the fitted population clearance is 5.00 L/h and its approximate uncertainty distribution is normal on the log scale, with standard deviation \(0.09531\). Take one hypothetical standard-normal draw \(z=1.00\). Hold \(V_{\mathrm{pop}}=25.0\) L fixed in this simplified example. The replicate's clearance draw is

\log(CL_{\mathrm{pop}}^{\mathrm{rep}}) =\log(5.00)+0.09531(1.00), \qquad CL_{\mathrm{pop}}^{\mathrm{rep}}=5.00\exp(0.09531)=5.50\ \text{L/h}.

For a 200 mg bolus and a time of 2.00 h, consider the structural concentration prediction with no individual effect in this arithmetic illustration and residual error omitted:

C^{\mathrm{rep}}(2.00) =\frac{200}{25.0}\exp\left(-\frac{5.50}{25.0}\times2.00\right) =8.00\exp(-0.440) =5.15229\ \text{mg/L}.

The new clearance applies to all individuals in this replicate. A new parameter vector can be drawn for the next replicate, but the draw is not repeated independently for each person. In a complete replicate, individual random effects and residual errors would also be generated unless the target conditions on them or the calculation is intentionally limited to this parameter-uncertainty illustration.

When to Simulate Doses and Observation Times

Doses and sampling times are part of the design, not automatically random variables. If the predictive question is whether the model reproduces concentrations under the actual study design, hold the observed dose histories and observation times fixed. This makes observed and replicated datasets comparable at the same design points. Simulate design quantities only when the predictive target includes their variation, and then use a stated design mechanism. For example, if dose amounts are assigned randomly by protocol, or observation times arise from a specified sampling process, those mechanisms can be included in the simulation.

Worked Example: A Random Observation Time

Suppose, hypothetically, a design specifies that one concentration is collected at a time drawn uniformly between 1.00 and 3.00 h after a 120 mg IV bolus. In this illustrative replicate, the time draw is \(t=2.00\) h. Use \(CL_{\mathrm{pop}}=6.00\) L/h, \(V_{\mathrm{pop}}=30.0\) L, and a zero individual-effect draw for a simple arithmetic illustration. The concentration prediction is

C(2.00)=\frac{120}{30.0}\exp\left(-\frac{6.00}{30.0}\times2.00\right) =4.00\exp(-0.400) =2.68128\ \text{mg/L}.

Let the log-scale residual standard deviation be \(0.150\), with a hypothetical standardized residual draw \(z=-1.00\). Then \(\epsilon=-0.150\), and the simulated observation is

Y^{\mathrm{rep}}=2.68128\exp(-0.150)=2.30780\ \text{mg/L}.

The time draw is a design draw, the zero individual effect is used only to keep this calculation transparent, and the residual draw creates observation-level variation. In an actual population-level replicate, individual effects would ordinarily be drawn from the model. If the target instead conditions on the observed time, use that time rather than drawing from the uniform design.

Diagnostics and Common Failure Modes

Before interpreting a predictive check, verify that the simulation actually follows the intended model. Useful checks include confirming parameter units and transformations, checking that all generated parameter values are admissible, and comparing the marginal distribution and correlation of simulated random effects with the specified \(\Omega\). For each simulated individual, verify that calculated concentrations respond sensibly to dose, time, clearance, and volume. Under the bolus model above, for example, concentration at time zero is \(D_i/V_i\), and concentration declines with time when clearance and volume are positive.

Wrong level for a draw. Drawing a new population parameter for every person inflates between-subject variation. Draw a population vector once per replicate, then draw individual effects separately for individuals.
Omitted residual error. Replicating only model predictions produces a smoother and less variable dataset than the observation model specifies. Add residual error when the predictive target concerns observed concentrations.
Incorrect error scale. A residual draw on the log scale must enter through the model's log-scale transformation, as in \(C\exp(\epsilon)\); it must not be added directly in concentration units.
Accidental design mismatch. Simulating new times or doses when the check is meant to condition on the observed design can make differences reflect design changes rather than model performance.
Broken parameter dependence. Drawing correlated population estimates or random effects independently can distort the joint distribution. Use the relevant covariance structure.
Confusing a conditional check with a population simulation. Re-drawing random effects for a target that is meant to condition on subject-specific information changes the question. State whether effects are newly drawn, held at estimates, or drawn from a conditional distribution.

A plug-in predictive simulation is often simpler and isolates the implications of the fitted model at its estimated parameter values. Incorporating parameter uncertainty can better represent uncertainty in future data predictions, but depends on a defensible uncertainty distribution and can broaden predictive variation. Conditioning on individual effects can answer a different, subject-specific question from generating new individuals from the population. None of these approaches is universally preferable: the simulation must match the predictive target defined for the check.

Common mistake: “I have simulated the model predictions, so I have generated replicate observations.” This is only true if the target is explicitly the distribution of predictions. For replicated concentration data under an observation model, the residual-error step is also required; and for a population-level target, between-subject variability must be represented as well.

Check Your Understanding

Answer each question by identifying the draw level or calculation required. Solutions follow.

In a population-level replicate, should a population clearance draw be shared by all individuals or drawn independently for each individual?
For a fixed-design predictive check, should observed sampling times ordinarily be held fixed or newly generated?
With a lognormal residual model \(Y=C\exp(\epsilon)\), \(C=4.00\) mg/L and \(\epsilon=0.10\), calculate the simulated observation to three significant figures.
Give one reason a plug-in simulation and a simulation with parameter uncertainty can produce different predictive variation.
If a replicate contains model-predicted concentrations but no residual error, what component of observed-data variation has been omitted?

Solutions. (1) Share one population clearance draw across individuals in that replicate; individual variation is generated through individual random effects. (2) Hold the observed times fixed when the target conditions on the observed design. (3) \(Y=4.00\exp(0.10)=4.42068\), or \(4.42\) mg/L to three significant figures. (4) The parameter-uncertainty simulation varies the population parameter vector across replicates, whereas the plug-in simulation fixes it at its estimate. (5) Observation-level residual variation specified by the residual-error model has been omitted.