Why Simulate a Fitted Mixture?
A fitted fast–slow clearance mixture can look numerically stable while still giving unreliable estimates. The fast class may be rare, the class-specific clearances may be difficult to distinguish at the available sampling times, or individual concentration profiles may not contain enough information for confident class assignment. Simulation checks ask a practical question: if the fitted model were a reasonable description of the population, would this study design let us recover the quantities we care about?
The central idea is to simulate many datasets from the fitted model, keeping the actual study design, and refit the same model to each dataset. The simulated class labels and parameter values are known, so estimates can be compared with their generating values. Class assignments can also be evaluated against the known simulated labels. These checks do not prove that the model is true; they reveal consequences of the fitted model and design.
Model and Simulation Assumptions
For a concrete illustration, suppose each individual belongs to a fast-clearance class \(F\) or a slow-clearance class \(S\). Let \(Z_i\) denote individual \(i\)'s class, with \(P(Z_i=F)=\pi_F\) and \(P(Z_i=S)=1-\pi_F\). A simple clearance model is \(CL_i=CL_{Z_i}\exp(\eta_i)\), where \(\eta_i\sim N(0,\omega^2)\). Thus, \(CL_F\) and \(CL_S\) are the class-specific typical clearances, and \(\omega^2\) describes within-class variability in log clearance.
For an illustrative one-compartment intravenous bolus model with fixed volume \(V\) and dose \(D\), the concentration predicted for an individual at time \(t\) is \(C_i(t)=\frac{D}{V}\exp[-CL_i t/V]\). To represent proportional lognormal residual error, let an observed concentration be \(Y_{ij}=C_i(t_{ij})\exp(\epsilon_{ij})\), where \(\epsilon_{ij}\sim N(0,\sigma^2)\). Other pharmacokinetic models and residual-error structures can be used; the simulation must use the model actually fitted.
The procedure assumes that the fitted model is the data-generating model for the purpose of the check, that its parameter values are suitable generating values, and that its random-effect and residual distributions are correctly specified. It also assumes that the design is represented correctly. In particular, retain the study's doses, sample times, covariates, and observation pattern for each individual. If the observed design includes missing samples, preserve that pattern for a conditional design check, or simulate missingness using an explicit missingness model when missingness itself is part of the question.
From Fitted Model to Replicate Datasets
Use the fitted parameter estimates as generating values. For every replicate dataset, draw each individual's class from the fitted class prevalence, draw their random effects from the fitted distribution, and generate observations at that individual's recorded design points. Then refit the same model, including the same class ordering and relevant constraints. Repeating the full fitting procedure matters: merely generating outcomes around fixed individual predictions would not test recovery of population class prevalence and clearances.
Record the fitted parameters and each individual's doses, covariates, sampling times, and observation pattern.
Draw class membership, random effects, and residual errors, then calculate observations at the recorded design points.
Apply the original model and fitting strategy to each replicate, noting convergence failures and boundary estimates as well as successful fits.
Compare fitted estimates with generating values and compare posterior class probabilities with the known simulated classes.
A basic measure of parameter recovery is the empirical bias across \(R\) successful replicate fits: \(\text{Bias}(\hat\theta)=R^{-1}\sum_{r=1}^{R}\hat\theta_r-\theta_0\), where \(\theta_0\) is the generating value. Also examine the spread of estimates, intervals, convergence rate, and whether estimates accumulate near parameter boundaries. Report how many fits were attempted and how many were included in each summary; a reassuring average among a small selected subset can hide systematic fitting failures.
Worked Examples
Worked Example: Simulating Classes at the Observed Sampling Time
Suppose a hypothetical study gives a 100 mg intravenous bolus, has \(V=20\) L, and samples every individual at \(t=2\) h. The fitted generating values are \(CL_F=5\) L/h, \(CL_S=2\) L/h, \(\pi_F=0.35\), and \(\sigma=0.15\). For this small illustration, set \(\eta_i=0\) for all individuals; a full simulation would draw \(\eta_i\) using the fitted \(\omega^2\).
The model-predicted concentrations are \(100/20=5\) mg/L before elimination, so for the fast class \(C_F(2)=5\exp[-(5)(2)/20]=5\exp(-0.5)=3.033\) mg/L. For the slow class, \(C_S(2)=5\exp[-(2)(2)/20]=5\exp(-0.2)=4.094\) mg/L.
For six hypothetical individuals, suppose simulated uniform draws for membership are \(0.12, 0.41, 0.83, 0.29, 0.67,\) and \(0.34\). Assign an individual to \(F\) when the draw is less than \(0.35\), and to \(S\) otherwise. The draws \(0.12, 0.29,\) and \(0.34\) give three fast-class individuals, so the simulated fast-class proportion is \(3/6=0.500\). This does not indicate an error: a small replicate can differ substantially from the generating prevalence of \(0.35\).
If the first individual's residual draw is \(\epsilon=0.10\), their simulated observation is \(3.033\exp(0.10)=3.352\) mg/L. If the second individual's residual draw is \(-0.05\), their slow-class observation is \(4.094\exp(-0.05)=3.894\) mg/L. These are simulated outcomes at the study's actual time point, not new sampling times.
The variability in simulated prevalence has a simple justification. For \(n\) independently assigned individuals, each fast-class indicator has expectation \(\pi_F\) and variance \(\pi_F(1-\pi_F)\). The mean of these indicators, the simulated proportion, therefore has expectation \(\pi_F\) and variance \(\pi_F(1-\pi_F)/n\). With \(\pi_F=0.35\) and \(n=6\), its standard deviation is \(\sqrt{0.35(0.65)/6}=0.195\). A single replicate is consequently a noisy check of prevalence; many replicates are needed to characterize its behavior.
Worked Example: Summarizing Parameter Recovery
Consider a hypothetical recovery exercise with 200 simulated datasets generated using \(\pi_F=0.35\), \(CL_F=5.00\) L/h, and \(CL_S=2.00\) L/h. Suppose 188 fits converge and the mean estimates among those converged fits are \(\bar{\hat\pi}_F=0.36\), \(\bar{\widehat{CL}}_F=5.12\) L/h, and \(\bar{\widehat{CL}}_S=1.96\) L/h.
The estimated biases are \(0.36-0.35=0.010\) for fast-class prevalence, \(5.12-5.00=0.12\) L/h for fast clearance, and \(1.96-2.00=-0.04\) L/h for slow clearance. Relative to the generating values, the clearance biases are \(0.12/5.00=2.4\%\) and \(-0.04/2.00=-2.0\%\), respectively. The convergence proportion is \(188/200=0.940\), or 94.0%.
These summaries suggest modest average bias among converged fits, but they do not by themselves establish good recovery. The 12 failed fits must be investigated, and the spread of estimates and interval coverage should be examined. For example, if failures occur mostly when the fast-class estimate is very small or the two estimates approach one another, the average of the remaining fits may conceal weak identification. The results are conditional on the selected generating parameters and design.
Checking Individual Class Uncertainty
For an observed individual, the posterior fast-class probability is the probability of membership in \(F\) conditional on their observations and the fitted model. It is proportional to the prior class probability multiplied by that class's likelihood for the individual's data. In an NLME model, the class likelihood ordinarily integrates over that individual's random effects. A single-point calculation with no random effect is a useful transparent illustration, but it is not a substitute for that integration in a full model.
Worked Example: A Posterior Class Probability From One Observation
Use the previous hypothetical study, with \(\pi_F=0.35\), class predictions \(3.033\) mg/L and \(4.094\) mg/L at two hours, and log-scale residual standard deviation \(\sigma=0.15\). Suppose one individual has an observed concentration of \(3.50\) mg/L. For this calculation, assume \(\eta=0\) and use only this observation.
Ignoring a factor common to both class likelihoods, the relative likelihood for a class prediction \(C\) is \(\exp\{-[\ln(y/C)]^2/(2\sigma^2)\}\). The fast-class log discrepancy is \(\ln(3.50/3.03265)\approx0.1433\), giving relative likelihood \(\exp[-(0.1433)^2/(2(0.15)^2)]\approx0.634\). The slow-class discrepancy is \(\ln(3.50/4.09365)\approx-0.1567\), giving relative likelihood \(\exp[-(-0.1567)^2/(2(0.15)^2)]\approx0.580\).
The posterior fast probability is therefore approximately \((0.35)(0.634)/\{(0.35)(0.634)+(0.65)(0.580)\}=0.2219/(0.2219+0.3770)=0.371\). The result is not a decisive assignment: under this simplified calculation, the fitted model gives the individual a 0.371 posterior probability of the fast class and a 0.629 probability of the slow class.
Across simulated individuals, compare posterior probabilities with the known simulated labels. Useful summaries include the distribution of probabilities among true fast and true slow individuals, the proportion of individuals receiving intermediate probabilities, and calibration: among individuals assigned similar probabilities, does the observed fraction of simulated fast labels roughly match those probabilities? A class probability near 0.5 is not a fitting failure by itself; it can correctly reflect limited information under the observed design.
Diagnostics, Alternatives, and Limitations
A simulation check is most informative when it examines more than a single mean estimate. Look for repeated nonconvergence, implausible or boundary estimates, highly skewed estimate distributions, strong dependence between class prevalence and class-specific clearance estimates, and poor recovery of class probabilities. Also check whether the generating parameter values are scientifically plausible and whether conclusions change when reasonable fitted values are used. If only some individuals have informative sampling times, summaries by design pattern can reveal which parts of the design contribute little information.
Simulation complements, rather than replaces, other checks. Standard errors and likelihood-based intervals summarize uncertainty under their approximations; simulation directly probes repeated behavior under specified generating assumptions. Visual diagnostics can reveal mismatches between observed and model-generated data, whereas a recovery exercise asks whether parameters could be estimated from data generated by the assumed model. Neither guarantees that the assumed mixture is appropriate for real data. This check is also distinct from deciding whether a second class should be included: here, the fitted mixture is treated as the generating model to examine its recovery properties.
- Do not change the design inadvertently. Simulating a new balanced schedule can make recovery look better than it would be for the observed, possibly sparse schedule.
- Do not interpret one replicate's prevalence as a parameter estimate. Class counts vary randomly around the generating prevalence, particularly in small samples.
- Do not omit failed fits. Report failures and inspect their patterns; summaries based only on successful fits can be misleading.
- Do not treat posterior probabilities as known labels. Uncertain probabilities are expected when individual data are weakly informative.
Check Your Understanding
Answer each question using the definitions and calculations above.
- If the observed study has \(n=40\) individuals and \(\pi_F=0.25\), what are the expected simulated fast-class proportion and its standard deviation?
- Why should a design-preserving simulation retain each individual's observed sample times?
- A recovery exercise has 100 attempted fits, 90 converged fits, and a mean fast-clearance estimate of \(4.80\) L/h when the generating value is \(5.00\) L/h. What are the convergence proportion and the estimated bias among converged fits?
- What does a posterior fast-class probability of 0.50 mean, and what does it not mean?
Solutions. (1) The expected proportion is \(\pi_F=0.25\). Its standard deviation is \(\sqrt{0.25(0.75)/40}=\sqrt{0.0046875}=0.0685\). (2) Sample times affect which concentration differences are observed and therefore the information available to distinguish classes; changing them would check a different design. (3) The convergence proportion is \(90/100=0.900\), or 90.0%. The bias among converged fits is \(4.80-5.00=-0.20\) L/h. (4) It means that, conditional on the fitted model and this individual's observations, the posterior probability of fast-class membership is 0.50. It does not mean the class is objectively half-fast, nor does it establish that the model is correct.