From Bootstrap Refits to Prediction Uncertainty
A pharmacometric prediction is often a function of several estimated parameters. Reporting only the prediction at the original fitted estimates can conceal how much it depends on the finite data used to fit the model. A bootstrap provides an empirical route to propagate that estimation uncertainty: resample the data according to the study design, refit the model, and evaluate the same prediction using each retained fit.
The resulting distribution is an approximation to the repeated-sampling uncertainty of the prediction under the fitted model and the chosen resampling scheme. It is not automatically a distribution of future patient outcomes. To include individual variability, residual error, or other sources of variation, those must be simulated separately in a way appropriate to the target. Here the focus is the uncertainty arising from parameter estimation.
A replicate is retained when it follows the prespecified resampling and fitting procedure and produces a usable estimate for the target. Retention should depend on documented procedural criteria, such as whether the optimizer converged and whether required estimates and predictions are finite—not on whether an estimate or prediction looks favorable or is near the original estimate. A valid estimate at a parameter boundary is not automatically a failed fit; it should be treated according to the model, estimation procedure, and prespecified rules.
Assumptions and a Practical Workflow
Bootstrap uncertainty is informative only when the resampling mechanism represents the data-generating design relevant to the target. For independent subjects, this generally means resampling subjects as units, carrying each selected subject’s records together. For clustered, repeated-measures, or otherwise dependent observations, the resampling unit must preserve the relevant dependence. Resampling individual observations as if they were independent when they are not can understate uncertainty or distort it.
The fitting procedure must also be applied consistently across replicates. Use the same model, estimation settings, parameter definitions, and target calculation. The bootstrap approximation further relies on the fitted model being sufficiently identifiable and on the resampling distribution of the estimates being a useful approximation to their sampling distribution. It can be unreliable with sparse data, weakly identified parameters, boundary behavior, or a high proportion of failed fits.
Define the prediction function, its units, and the covariate or dosing scenario. Hold that target fixed across replicates unless the estimand explicitly requires it to vary.
Generate bootstrap data sets using the independent sampling unit or a design-appropriate scheme. Record the random seed and resampling settings.
Store each replicate’s parameter vector, convergence information, warnings, and any failure reason. Do not silently replace failed fits or discard extreme estimates.
For each retained parameter vector, calculate the prediction using the same target definition and parameter transformations as for the original fit.
Calculate the chosen empirical quantiles, report the number and fraction of successful fits, and check whether failures or Monte Carlo variability undermine interpretation.
Why the Empirical Interval Works—and What It Means
Let \(B\) be the number of retained bootstrap replicates. The empirical distribution of their predictions assigns probability \(1/B\) to each calculated value, counting repeated values with their multiplicity. Its empirical cumulative distribution function is
where \(I(\cdot)\) equals 1 when the condition is true and 0 otherwise. For a central \(100(1-\alpha)\%\) percentile interval, take the empirical quantiles at \(\alpha/2\) and \(1-\alpha/2\). Quantile algorithms differ in how they interpolate between ordered values, so the algorithm should be recorded; with a small \(B\), the choice and the limited resolution can matter.
The justification is a resampling approximation. If the bootstrap distribution of the refitted estimator approximates its repeated-sampling distribution, applying the same function \(g\) to every bootstrap estimate approximates the distribution of the corresponding prediction estimator. Taking its lower and upper quantiles then gives a percentile interval. This argument depends on the quality of the bootstrap approximation; the interval is not guaranteed to have its nominal coverage for every model or data set.
Worked Example: Retaining Replicates After Fit Failures
Suppose a hypothetical analysis requests 10 subject-level bootstrap refits. Eight converge and return finite estimates of clearance, in L/h: 4.6, 4.9, 5.0, 5.1, 5.2, 5.4, 5.7, and 6.1. Two runs stop without a usable estimate: one reaches the iteration limit and one returns a non-finite objective value. Assume the recorded diagnostics confirm that these are computational failures, not valid boundary solutions.
Retention and reporting. Retain the eight usable fits under the stated computational rule, while retaining the failed-fit records in the audit output. The success fraction is \(8/10=0.800\), or 80.0%. The mean of the eight retained clearance estimates is
The mean is only a simple summary here; it is not a substitute for examining the distribution or the fit diagnostics. Nor does 80.0% success establish that the retained estimates represent all bootstrap samples. With only 10 attempted replicates, the success fraction itself is imprecise, and two failures warrant investigation before relying on an interval. The analyst should inspect warnings and data sets, try justified numerical remedies consistently, and report how conclusions depend on the failed replicates. They should not choose to rerun only estimates that appear extreme or inconvenient.
If failures are rare and demonstrably computational, an interval based on the retained estimates can be presented with the failure rate and a limitation statement. If failures are frequent or systematically associated with particular resamples or parameter regions, the successful fits may be a selected subset. In that case, their quantiles need not represent the intended bootstrap distribution, so reporting a nominal interval without qualification is not justified.
Worked Example: Calculating a Prediction Interval from Retained Draws
Consider hypothetical retained bootstrap clearance estimates \(CL\) of 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, and 7.5 L/h. For a fixed one-compartment intravenous bolus scenario with dose \(D=500\) mg, volume \(V=50\) L, and time \(t=8\) h, the concentration prediction is
The units are mg/L because \(D/V=500\text{ mg}/50\text{ L}=10\text{ mg/L}\), while \(CLt/V\) is dimensionless. For the replicate with \(CL=5.0\) L/h,
Applying the same calculation to all ten estimates gives these predictions, rounded to three decimal places:
| CL (L/h) | Prediction C(8) (mg/L) |
|---|---|
| 3.0 | 6.188 |
| 3.5 | 5.712 |
| 4.0 | 5.273 |
| 4.5 | 4.868 |
| 5.0 | 4.493 |
| 5.5 | 4.148 |
| 6.0 | 3.829 |
| 6.5 | 3.535 |
| 7.0 | 3.263 |
| 7.5 | 3.012 |
For illustration, use the nearest-rank empirical quantile rule: the \(q\) quantile is the ordered observation at rank \(\lceil Bq\rceil\). There are \(B=10\) predictions. Ordered from smallest to largest, they are 3.012, 3.263, 3.535, 3.829, 4.148, 4.493, 4.868, 5.273, 5.712, and 6.188 mg/L. The 10th percentile is at rank \(\lceil10(0.10)\rceil=1\), and the 90th percentile is at rank \(\lceil10(0.90)\rceil=9\). Thus, the illustrative 80% percentile interval is \([3.012,\ 5.712]\) mg/L. Its coarse endpoints reflect the small number of retained draws; it is not a recommendation to use only 10 refits for a final analysis.
The median under the usual even-sample convention is the average of the fifth and sixth ordered predictions: \((4.148+4.493)/2=4.321\) mg/L. This is the median bootstrap prediction, not necessarily the prediction evaluated at the mean parameter estimates. Since concentration decreases as clearance increases in this model, sorting the predictions reverses the clearance ordering. Quantiles should therefore be calculated from the prediction draws themselves, rather than assuming the parameter and prediction endpoints have the same order.
Worked Example: Propagating Paired Parameter Estimates
A bootstrap refit usually returns a vector of estimates, not independent estimates for each parameter. Consider four hypothetical retained pairs \((CL,V)\), in L/h and L, respectively: (4, 40), (5, 40), (5, 60), and (7, 60). For the same 500 mg bolus and \(t=8\) h, calculate \(C(8)\) from each pair, preserving the values from the same refit.
For \((CL,V)=(4,40)\), \(D/V=12.5\) mg/L and \(CLt/V=4(8)/40=0.8\), so \(C(8)=12.5\exp(-0.8)=5.617\) mg/L. For \((5,40)\), the exponent is \(-5(8)/40=-1\), giving \(12.5\exp(-1)=4.598\) mg/L. For \((5,60)\), the prefactor is \(500/60=8.333\) mg/L and the exponent is \(-5(8)/60=-0.6667\), giving \(8.333\exp(-0.6667)=4.278\) mg/L. For \((7,60)\), the exponent is \(-7(8)/60=-0.9333\), giving \(8.333\exp(-0.9333)=3.277\) mg/L.
The four predictions are therefore 5.617, 4.598, 4.278, and 3.277 mg/L. Their mean is \((5.617+4.598+4.278+3.277)/4=4.443\) mg/L, subject to rounding. With only four draws, the minimum-to-maximum range is 3.277 to 5.617 mg/L; it should not be presented as a stable nominal percentile interval.
Each pair represents the joint estimate from one bootstrap refit. Recombining clearance from one replicate with volume from another would create parameter vectors that the bootstrap did not produce and could change the prediction distribution. The appropriate unit of propagation is the retained parameter vector, with its within-replicate relationships intact.
Diagnostics, Alternatives, and Limitations
A reproducible result should make it possible to reconstruct the route from each resample to each prediction. Store the replicate identifier, resampling information, fit status and warnings, retained parameter vector, and prediction. Report the attempted and retained counts, the rule for declaring a fit usable, the quantile method, the target scenario, and the number of digits used in the calculation. Do not round parameter estimates before evaluating a nonlinear prediction; round for display after calculating the predictions and interval.
Useful diagnostics include the distribution of each parameter and prediction, the locations of failed fits, and the stability of interval endpoints as more replicates are generated. If endpoints change materially across independent runs or when the number of replicates increases, Monte Carlo error may be consequential. If optimizer warnings, boundary estimates, or failures cluster in particular regions, that pattern may indicate weak identification or a model-fitting problem rather than a harmless numerical nuisance. More bootstrap replicates can reduce Monte Carlo noise, but cannot repair a resampling scheme that fails to represent the design or a model that is poorly identified.
A common alternative is to draw parameters from a multivariate normal approximation based on the estimated covariance matrix, then evaluate the same prediction function. This can be efficient when the estimator is approximately normal and the covariance estimate is reliable, but may poorly represent skewness, boundaries, or irregular sampling distributions. Bootstrap refits use the empirical resampling distribution instead, at the cost of repeated model fits and possible convergence failures. A Bayesian simulation instead propagates posterior parameter draws and has a different inferential interpretation; it is not interchangeable with a bootstrap sampling distribution.
The percentile interval is also not the only bootstrap interval construction. Other constructions adjust for bias or the shape of the sampling distribution, but require additional assumptions and calculations. For the percentile method taught here, the main safeguards are to define the estimand and resampling scheme clearly, avoid value-based selection of replicates, preserve joint parameter vectors, and report limitations imposed by fit failures and finite \(B\).
Common Pitfalls
- Dropping extreme estimates because they seem implausible. This conditions the empirical distribution on the analyst’s preferred range and can narrow or shift the interval. Investigate whether an estimate reflects a valid fit or a documented computational failure; do not silently filter by its numerical value.
- Calling every boundary estimate a failed fit. A boundary result can be a valid outcome of the specified estimation procedure. Assess and report it using prespecified criteria rather than treating it as a numerical failure solely because it is inconvenient.
- Ignoring failures when computing quantiles. Quantiles of successful fits alone may be selection-biased if failure depends on the resampled data or estimate. Report failures and qualify or withhold the interval when success is inadequate.
- Calling parameter-only uncertainty a full outcome prediction interval. The interval above does not include residual or individual variability. State exactly which sources of uncertainty the simulation includes.
- Breaking up parameter vectors. Independently resampling parameter columns discards the dependence represented by the paired bootstrap estimates and can produce a different prediction distribution.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A bootstrap analysis has 200 attempted refits, of which 184 meet a prespecified convergence and finite-output rule. What success fraction should be reported, and what else should accompany it?
- For a fixed target, retained prediction draws in mg/L are 2.0, 2.5, 3.0, 3.5, and 4.0. Under the nearest-rank rule, what are the 20th and 80th percentile values?
- Why should a prediction based on a bootstrap replicate use the full parameter vector from that same replicate?
- Does an interval formed only by applying a deterministic concentration model to parameter draws include residual error or between-subject variability? Explain.
Solutions. (1) The success fraction is \(184/200=0.920\), or 92.0%. Report the failure count and reasons, the retention rule, and diagnostics sufficient to assess whether failures may select the retained fits. (2) The ordered values are already shown; ranks are \(\lceil5(0.20)\rceil=1\) and \(\lceil5(0.80)\rceil=4\), so the endpoints are 2.0 and 3.5 mg/L. (3) The paired estimates represent the joint result of one refit. Mixing parameter values across refits changes that joint distribution and may alter the prediction uncertainty. (4) No. It includes parameter-estimation uncertainty propagated through the deterministic model, but residual and between-subject variability require additional simulation if the target includes them.