Why Compare Individual Sampling Schedules?
An individual pharmacokinetic (PK) schedule specifies when concentrations are measured for one person. With sparse sampling, each measurement has to contribute useful information: measurements clustered in time may describe concentration levels but provide little information about the rate of elimination. Simulation lets us compare candidate schedules under a specified PK model and assess the resulting estimation bias, precision, and confidence-interval coverage before treating one schedule as preferable.
The comparison below focuses on clearance in a one-compartment model after an intravenous bolus dose. It is an individual-level comparison: each simulated person receives one schedule, and the target is that person's clearance. The calculations and simulation results are hypothetical. They illustrate a comparison under stated assumptions; they are not results from a clinical study.
Model, Assumptions, and Information in the Sampling Times
Suppose concentration follows a one-compartment exponential decline. On the log scale, write the observation at time \(t_j\) as
Here \(Y_j=\log(C_j)\), \(C_j\) is concentration, \(a\) is an intercept, \(CL\) is clearance in L/h, \(V\) is volume in L, and \(m\) is the number of samples for that individual. The intercept is estimated rather than treated as known; it can absorb an unknown concentration scale. Assume \(V\) is known, and the errors are independent with \(\epsilon_j\sim N(0,\sigma^2)\), where \(\sigma\) is the standard deviation of log-scale measurement error. The times are fixed by the design, and the linear model is correctly specified.
These assumptions matter. A systematic time trend in residuals, correlated errors, concentration-dependent error not represented by the model, or a misspecified elimination profile can make the stated precision calculation unreliable. The derivation also assumes that all planned concentrations are observed and usable. Missing samples and measurements below quantification need to be represented in the simulation if they are expected in practice.
The slope of the regression of \(Y_j\) on \(t_j\) estimates \(-CL/V\). Let \(S_{tt}=\sum_{j=1}^{m}(t_j-\bar{t})^2\), where \(\bar{t}\) is the mean sampling time. Since the error variance is \(\sigma^2\), the variance of the estimated slope is \(\sigma^2/S_{tt}\). Multiplying the slope estimate by \(-V\) gives the clearance estimate, so its standard error is \(V\sigma/\sqrt{S_{tt}}\). Thus it is the spread of the sampling times around their mean, not simply the latest sampling time, that governs precision in this model.
Worked Example: Comparing the Information in Two Schedules
Consider a hypothetical individual with known \(V=40\) L and log-scale error standard deviation \(\sigma=0.10\). Schedule A samples at 0.5, 1.5, 4.5, and 5.5373 h. Schedule B samples at 0.5, 2, 5, and 8 h. Both schedules have four samples, so this comparison holds the number of samples fixed.
For Schedule A, the sampling-time mean is approximately 3.0093 h. The centered sum of squares is approximately \(17.1875\ \text{h}^2\). Therefore,
For Schedule B, \(\bar{t}=(0.5+2+5+8)/4=3.875\) h. Its centered sum of squares is
Thus \(SE_B=4/\sqrt{33.1875}\approx0.6943\ \text{L/h}\). In this model, B is more precise because its sampling times are more spread out. This is a conditional model-based comparison, not proof that B will perform better under every real data-generating process.
For a known error standard deviation, a 95% normal interval for clearance has half-width \(1.959964\,SE\). Schedule A's half-width is \(1.959964\{4/\sqrt{17.1875}\}\approx1.8910\ \text{L/h}\), giving a total width of about \(3.7821\ \text{L/h}\). Schedule B's half-width is approximately \(1.3609\ \text{L/h}\), giving a width of about \(2.7218\ \text{L/h}\). A narrower interval indicates greater precision under this model; it does not establish that the interval procedure has adequate coverage in a more complex setting.
Simulate the Comparison and Summarize Performance
A useful comparison specifies the target parameter, simulation conditions, estimator, interval procedure, and performance measures before looking at results. For each replicate, generate concentrations under the same true clearance and volume for both schedules, add measurement errors under the assumed error model, estimate clearance, and construct the prespecified interval. Repeat this process many times. If a replicate fails to produce an estimate, record and report that failure rather than silently dropping it.
For \(R\) replicates, let \(\widehat{CL}_r\) be the estimate in replicate \(r\), and let \(CL_{\mathrm{true}}\) be the simulated true clearance. The estimated bias is the average signed estimation error, and the root mean squared error (RMSE) combines bias and variability. Coverage is the proportion of intervals containing the simulated true value. These summaries answer different questions: low bias does not guarantee low variability, and a narrow interval does not guarantee correct coverage.
Here \([L_r,U_r]\) is the interval from replicate \(r\), and \(I\{\cdot\}\) equals 1 when the condition inside braces is true and 0 otherwise. Coverage is an estimate of the procedure's repeated-sampling probability of containing the true clearance under the simulated conditions. It is not the probability that a particular interval contains the parameter.
Worked Example: Reading Bias, RMSE, and Coverage Together
Suppose a hypothetical simulation uses \(CL_{\mathrm{true}}=4.00\) L/h and \(R=1000\) replicates per schedule. For Schedule A, the average estimate is 4.012 L/h, the empirical standard deviation of estimates is 0.9700 L/h, and 947 of 1000 nominal 95% intervals contain the true value. For Schedule B, the average estimate is 3.991 L/h, the empirical standard deviation is 0.7000 L/h, and 951 intervals contain the true value.
The estimated biases are \(4.012-4.00=0.012\) L/h for A and \(3.991-4.00=-0.009\) L/h for B. Their signs indicate slight average overestimation for A and underestimation for B in these simulated replicates. The bias is small relative to the empirical spread in each case.
The RMSE can also be checked from the identity \(\mathrm{RMSE}^2=\mathrm{Bias}^2+\mathrm{Var}(\widehat{CL})\), using the replicate distribution. Using the reported empirical standard deviations as the spread estimates gives approximately
The identity uses the variance, not the standard deviation, and the displayed approximation uses the reported rounded summaries. The observed coverage rates are \(947/1000=94.7\%\) for A and \(951/1000=95.1\%\) for B. Both are close to the nominal 95% in this run. The small difference in rates alone is not persuasive evidence that one interval procedure has better coverage; finite simulation adds uncertainty to each estimated rate.
Paired Comparisons and Their Uncertainty
When the same simulated individuals and underlying random conditions are used for both schedules, compare replicate-level outcomes rather than treating the two sets of results as unrelated. This is the paired-comparison principle. For example, let \(D_r\) be the difference between Schedule A's and Schedule B's absolute estimation errors for the same simulated replicate. The mean of \(D_r\) estimates the average difference in absolute error; a positive value favors B on that metric. The standard error of the average paired difference is the sample standard deviation of the \(D_r\)'s divided by \(\sqrt{R}\).
This comparison requires correctly matched replicates: each difference must use outcomes generated under the same replicate conditions. Common random numbers can make the paired differences less variable when they induce positive correlation between the two outcomes. They do not guarantee a reduction in uncertainty, so the observed variability of the paired differences should be checked.
Worked Example: Uncertainty in the Paired Error Difference
Suppose the 1000 paired replicates in the hypothetical comparison produce a mean difference in absolute clearance error of 0.190 L/h, defined as error under A minus error under B, with a sample standard deviation of paired differences of 0.310 L/h. The estimated standard error is
Using a large-sample 95% interval, \(\bar{D}\pm1.96SE(\bar{D})\), gives \(0.190\pm1.96(0.00980)\), or approximately \([0.1708,0.2092]\) L/h. Under this definition, the interval suggests lower average absolute error for B in the simulated conditions. It quantifies Monte Carlo uncertainty in the average paired difference; it does not measure uncertainty about how well the simulation model represents clinical practice.
Diagnostics, Alternatives, and Limitations
A comparison should be checked before a numerical ranking is used. First, verify that each schedule is implemented at the intended times and that the same true parameter values and estimation rules are used. Inspect convergence and estimation failures, the distribution of estimates across replicates, and interval widths and coverage. For the log-linear model, residual plots against time can reveal curvature or changing variance. In a nonlinear estimator, boundary estimates, skewed estimation errors, or many failed fits can make a mean and standard deviation incomplete summaries.
The analytic standard error derived above is a useful benchmark when its assumptions hold. Simulation is more flexible: it can represent nonlinear estimation, censoring, missing samples, and more realistic error structures, if those features are implemented. Conversely, simulation can give precise answers to the wrong model. A reasonable analysis therefore checks performance under plausible deviations from the primary assumptions and reports which conditions were simulated.
Different summaries may lead to different design choices. Bias emphasizes systematic error, RMSE combines bias and spread, and coverage evaluates interval calibration. Precision alone may favor an informative schedule that is operationally difficult or yields unusable late samples. A fixed number of samples makes the comparison easier to interpret, but schedules with different sample counts may also be compared when burden or cost is part of the design question. State what is held constant and why.
A common mistake is to choose the schedule with the latest sample because it seems to capture elimination best. In the stated model, precision depends on \(S_{tt}\), the total spread around the mean, so one late sample does not automatically outweigh several measurements concentrated early. The correction is to calculate or simulate performance for the complete schedule, then assess whether the model assumptions and operational constraints support the resulting comparison.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Why does the standard error formula use centered sampling times when the intercept is estimated?
- Under the model in this tutorial, which schedule has greater model-based precision: one with \(S_{tt}=12\ \text{h}^2\) or one with \(S_{tt}=30\ \text{h}^2\)? Explain.
- A simulation reports bias of 0.20 L/h and standard deviation of 0.50 L/h. What is the corresponding RMSE?
- What does a simulated coverage rate of 93% mean for a nominal 95% interval procedure?
Solutions. 1. Estimating the intercept removes the component shared by all observations; slope precision depends on deviations \(t_j-\bar{t}\), giving \(S_{tt}\). 2. The schedule with \(S_{tt}=30\ \text{h}^2\) has greater model-based precision because \(SE=V\sigma/\sqrt{S_{tt}}\) is smaller when \(S_{tt}\) is larger, all else equal. 3. \(\mathrm{RMSE}=\sqrt{0.20^2+0.50^2}=\sqrt{0.29}\approx0.5385\) L/h. 4. In the simulated conditions, about 93% of repeated intervals contained the true parameter. This is below the nominal 95%; it is an estimated coverage rate subject to Monte Carlo uncertainty, not a probability statement about any one interval.