Why Simulate BQL Handling Methods?
When concentrations fall below the lower limit of quantification (LLOQ), the value is not known precisely; the observation indicates that the concentration is below a threshold. The earlier material on left censoring and censored likelihoods establishes how this information can enter an analysis. A simulation study asks a different question: under a specified population PK model and censoring process, how well do alternative analysis methods recover known parameters?
A simulation makes the truth available because the investigator specifies the data-generating model and its parameter values. Repeating the process reveals how estimates vary across hypothetical data sets, how often uncertainty intervals include the truth, and how often estimation fails to converge. The conclusions are conditional on the simulated model and design; they are not universal rankings of methods.
For a fair comparison, define the data-generating process before fitting any method. For example, generate subject-specific random effects, generate concentrations with a specified structural and residual-error model, and record a value as BQL when its generated concentration is below the applicable LLOQ. The generated model and parameters are the simulation truth. Changing the LLOQ across scenarios can produce different censoring levels while leaving the true PK parameters unchanged.
In each replicate, analyze the same generated study using the methods of interest. A censored-likelihood method uses the probability of falling below the LLOQ for a BQL record. A missing-record approach omits BQL records from the concentration likelihood. A substitution approach replaces each BQL value with a chosen constant, such as LLOQ/2, and analyzes that substituted value as if it were observed. These choices encode different statistical treatments of the same information; none changes how the simulated data were generated.
Designing a Controlled Comparison
A useful design varies censoring level deliberately while keeping other conditions comparable. For instance, define low, moderate, and high censoring scenarios by selecting LLOQs that yield approximately 10%, 40%, and 70% BQL observations. These are targets, not guaranteed outcomes: the realized proportion varies among replicates. Report the realized censoring distribution as well as the target scenario.
The simulation should specify the population parameters to evaluate, such as clearance, volume, or a random-effect variance. It should also specify the number of subjects and samples, sampling times, random-effect and residual-error distributions, and the rule for generating BQL records. All methods should receive the same generated replicate. This paired design makes method differences easier to interpret because a particularly informative or difficult simulated data set affects all methods in that replicate.
The number of replicates, \(B\), must be large enough to estimate performance summaries with useful Monte Carlo precision. More replicates improve the precision of those summaries, but do not repair an unrealistic data-generating model. Record each replicate’s estimates, interval limits if available, convergence status, and realized BQL proportion. Do not silently discard failed fits.
Worked Example: Setting Up Censoring Scenarios
Suppose each hypothetical study has 100 subjects with six planned concentrations per subject, for 600 planned observations. The simulation uses the same population PK parameters and sampling schedule in three scenarios, with LLOQs chosen to target 10%, 40%, and 70% BQL observations.
If a replicate in the moderate-censoring scenario has 600 planned records and 246 are BQL, its realized BQL proportion is \(246/600=0.410\), or 41.0%. The target was 40%; the small difference is ordinary replicate-to-replicate variation, not a reason to relabel the scenario.
Analyze this same replicate by censored likelihood, missing-record deletion, and substitution. If the substitution rule is LLOQ/2, apply it to all 246 BQL records in that replicate. The resulting substituted data set is an analysis input, not a claim that the unobserved concentrations actually equaled LLOQ/2.
Performance Measures Across Replicates
Let \(\theta_0\) be a true scalar parameter and let \(\hat{\theta}_b\) be its estimate from replicate \(b\), for \(b=1,\ldots,B\). For summaries of estimates, first make clear which replicates are included. A common choice is to calculate estimate-based summaries among converged fits and report the convergence rate separately. Because convergence may depend on the data, these summaries describe estimates conditional on convergence and should not be presented without the failure rate.
Bias has the same units as the parameter. Empirical standard deviation and RMSE also have parameter units. For a parameter such as clearance measured in L/h, all three are in L/h. Relative bias may be reported as \(100(\bar{\hat{\theta}}-\theta_0)/\theta_0\) percent when \(\theta_0\ne 0\), but it does not replace reporting the parameter-scale result.
Precision here refers to the empirical spread of estimates across simulated data sets, not the width of one fitted interval. Coverage assesses whether the interval procedure achieves its nominal long-run inclusion frequency under the specified simulation. For example, 95% coverage means that, over repeated studies generated and analyzed under the stated conditions, approximately 95% of the intervals contain the fixed truth. It does not mean there is a 95% probability that the particular fixed parameter lies in a calculated interval.
Worked Example: Bias, Spread, and RMSE
In four hypothetical converged fits, estimated clearance values are 9, 10, 11, and 12 L/h. The true clearance is 10 L/h. The mean estimate is \((9+10+11+12)/4=10.5\) L/h, so empirical bias is \(10.5-10=0.500\) L/h.
The deviations from the mean are \(-1.5,-0.5,0.5,1.5\) L/h. Their squared sum is \(2.25+0.25+0.25+2.25=5.00\), giving sample standard deviation \(s=\sqrt{5/(4-1)}=1.291\) L/h.
The squared errors from the truth are \(1,0,1,4\), so RMSE is \(\sqrt{(1+0+1+4)/4}=\sqrt{1.5}=1.225\) L/h. The nonzero bias contributes to RMSE; the two summaries answer related but distinct questions.
Coverage, Convergence, and Monte Carlo Precision
For a method with fitted intervals \([L_b,U_b]\), define an indicator \(I_b=1\) when \(L_b\leq\theta_0\leq U_b\), and \(I_b=0\) otherwise. If all \(B\) replicates yield intervals, estimated coverage is \(C=B^{-1}\sum_{b=1}^{B} I_b\). Its approximate Monte Carlo standard error is \(\sqrt{C(1-C)/B}\), reflecting simulation noise in the estimated coverage.
For example, if 1,000 replicate intervals produce estimated coverage \(C=0.940\), the Monte Carlo standard error is approximately \(\sqrt{0.940(0.060)/1000}=0.00751\), or 0.751 percentage points. A small discrepancy from nominal 95% coverage may therefore be partly due to finite simulation size. Increasing \(B\) reduces this Monte Carlo error; it does not alter the underlying method or data-generating process.
Convergence requires an operational definition set before examining the results. It may include a successful optimizer termination, finite parameter estimates, and satisfaction of prespecified estimation diagnostics. Report the number and proportion of successful fits for every method and scenario. If some fits fail, state whether coverage is calculated among converged fits, and show that denominator. A convergence rate of 90% and coverage of 95% among converged fits do not mean that 95% of all attempted fits produced a covering interval.
The relationship between bias, spread, and RMSE can be justified directly. Define the empirical bias \(d=\bar{\hat{\theta}}-\theta_0\). Write each error as \(\hat{\theta}_b-\theta_0=(\hat{\theta}_b-\bar{\hat{\theta}})+d\), square and average over replicates. The cross-product sums to zero because \(\sum_b(\hat{\theta}_b-\bar{\hat{\theta}})=0\). Thus:
Here \(s^2\) is the usual sample variance with denominator \(B-1\). This identity shows why a method can have modest empirical spread but substantial RMSE when its estimates are systematically displaced from the truth.
Worked Example: Reporting Coverage With Failed Fits
Consider 100 planned fits for one method and scenario. Ninety-two converge. Among those 92 fits, 82 intervals contain the true parameter. The convergence rate is \(92/100=0.920\), or 92.0%. Conditional empirical coverage among converged fits is \(82/92=0.891\), or 89.1%.
The study should report both values, with the coverage denominator stated. It should not report “82% coverage” without explaining that this divides covered intervals by all 100 planned fits, including eight failures that produced no interval. Nor should it omit the eight failures and imply the method always converged. The distinction matters because the subset of converged replicates may not represent all simulated replicates.
With only 92 intervals, the approximate Monte Carlo standard error for the conditional coverage is \(\sqrt{0.891(1-0.891)/92}=0.0325\), or 3.25 percentage points. This relatively large value cautions against treating 89.1% as a precise estimate of the method’s long-run coverage.
Comparing the Methods and Diagnosing Problems
| Approach | What is done with a BQL record? | Simulation results to examine |
|---|---|---|
| Censored likelihood | Contributes the probability of a concentration below the LLOQ | Bias, spread, coverage, and convergence as censoring increases |
| Missing-record deletion | Is omitted from the concentration likelihood | Whether loss of low-concentration information changes estimates or interval performance |
| Substitution | Is replaced by a specified value and treated as observed | Dependence of results on the substitution rule and censoring level |
A censored likelihood can use the information that a concentration lies below a threshold, provided its probability model is appropriate and computation succeeds. Deletion discards that threshold information; substitution treats an invented value as if measured. A simulation can quantify consequences for the chosen design, but it cannot establish that one method is best in every PK model or assay setting.
Inspect results by parameter and censoring scenario, not just through one overall average. Increasing bias, deteriorating coverage, wider empirical spread, or falling convergence rates at higher censoring are warning signs that should be reported. A pattern in which estimates look stable only among a shrinking set of successful fits is especially important to flag. Also verify that the realized BQL proportions match the intended scenario broadly and that all methods were applied to the same replicate data.
A failure to recover the simulation truth does not by itself diagnose its cause. It may reflect information loss from censoring, a method’s treatment of BQL values, finite-sample variability, or numerical difficulty. Conclusions are limited to the specified data-generating assumptions, design, parameter values, interval procedure, and convergence criteria. This comparison is not a substitute for examining whether the simulation assumptions are credible for the application.
Common Pitfall
Correction: Small empirical spread means estimates varied less across the included replicates; it does not establish closeness to the truth. Compare bias and RMSE as well, assess interval coverage, and report convergence. A tightly clustered set of systematically biased estimates can be imprecise in the sense of error from truth despite having low empirical spread.
Another frequent error is treating a simulated BQL value substituted with LLOQ/2 as the true concentration. The value is an analysis convention, not information supplied by the observation. The simulation’s true generated concentration remains known to the investigator and can be used to define performance, but it is not available to the method analyzing that BQL record.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A simulation has 500 planned concentration records, of which 180 are BQL in one replicate. What is the realized censoring proportion?
- For a parameter with true value 20 mg/h, the mean estimate across converged replicates is 18 mg/h. What is the empirical bias, and what does its sign indicate?
- Why should a simulation report convergence rate separately from coverage among converged fits?
- A method has 95% empirical coverage from 200 intervals. Calculate the approximate Monte Carlo standard error of that coverage proportion.
- What information does a censored-likelihood analysis retain for a BQL record that substitution does not represent in the same way?
Solutions. (1) The realized proportion is \(180/500=0.360\), or 36.0%. (2) Bias is \(18-20=-2\) mg/h, indicating that the mean estimate is 2 mg/h below the truth. (3) Coverage among successful fits is conditional on convergence, while the convergence rate describes how often an interval-producing fit was obtained; failures may be selective. (4) The approximate standard error is \(\sqrt{0.95(0.05)/200}=0.0154\), or 1.54 percentage points. (5) The censored likelihood uses the probability that the concentration falls below the LLOQ; substitution inserts a chosen value and treats it as observed, rather than directly representing the below-threshold event.