Bootstrap Uncertainty for Derived PK Quantities

Why Derived PK Quantities Need Their Own Uncertainty Analysis

Population PK models are usually parameterized in quantities such as clearance, volumes, absorption rates, and residual variability. Clinical interpretation, however, often depends on quantities calculated from those parameters: exposure, half-life, or a volume-clearance ratio. A point estimate can be transformed directly, but its uncertainty cannot generally be obtained by transforming a standard error or interval in the same way.

The bootstrap provides a direct route: calculate the derived quantity separately for every successful bootstrap fit, then summarize the resulting distribution. The essential detail is to use the joint parameter estimates from each replicate. This retains the parameter dependence represented by the bootstrap, including dependence that may strongly affect a ratio or other nonlinear function.

Definition: Let \(\hat{\boldsymbol{\theta}}^{(b)}\) be the vector of parameter estimates from successful bootstrap replicate \(b\), for \(b=1,\ldots,B\). For a derived quantity \(g(\boldsymbol{\theta})\), calculate \(g^{(b)}=g(\hat{\boldsymbol{\theta}}^{(b)})\) within each replicate. The empirical distribution of \(g^{(1)},\ldots,g^{(B)}\) describes bootstrap uncertainty in that derived quantity.

The bootstrap replicates must come from a resampling scheme appropriate to the study design and must represent refitting the model when the aim is to quantify estimation uncertainty. The resampling unit and any stratification or clustering should be handled as specified by that design. The fitted model must also provide the parameters needed for the derived quantity, and the formula must be defined over the replicate estimates. For example, a ratio is not defined when its denominator is zero.

Transform the Replicates, Not Just the Summary

For each successful fit, apply the PK formula to that fit's parameter vector. The bootstrap standard error of the derived quantity is the sample standard deviation of the \(B\) transformed values. A bootstrap interval can likewise be formed from their empirical distribution, using an interval method appropriate to the analysis. The percentile interval, for example, uses the relevant empirical quantiles of the transformed values.

\widehat{SE}_{\mathrm{boot}}(g) = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B} \left(g^{(b)}-\bar{g}_{\mathrm{boot}}\right)^2}, \qquad \bar{g}_{\mathrm{boot}}=\frac{1}{B}\sum_{b=1}^{B}g^{(b)}.

This is the sample standard deviation across bootstrap replicates, not the standard deviation of individual observations within the original data. Each \(g^{(b)}\) is a model-derived value, and the resulting uncertainty concerns the estimated quantity under the chosen bootstrap procedure.

Worked Example: Bootstrap Uncertainty in Half-Life

Suppose six successful hypothetical bootstrap fits estimate the elimination rate constant \(k_e\). In each fit, half-life is calculated as \(t_{1/2}^{(b)}=\ln(2)/k_e^{(b)}\). The resulting six half-lives, in hours, are 5.0285, 5.4171, 5.8057, 6.1943, 6.5829, and 6.9715. These values are illustrative; six replicates are too few for a stable inferential interval.

Their mean is \(6.000\) h. Subtracting that mean gives deviations of approximately \(-0.9715\), \(-0.5829\), \(-0.1943\), \(0.1943\), \(0.5829\), and \(0.9715\) h. The bootstrap standard error is the sample standard deviation of these transformed values:

\widehat{SE}_{\mathrm{boot}}(t_{1/2}) = \sqrt{\frac{ 2(0.9715^2+0.5829^2+0.1943^2) }{6-1}} = \sqrt{\frac{2.6427}{5}} \approx 0.727\ \mathrm{h}.

Thus, in this illustrative set, the bootstrap standard error of half-life is about \(0.727\) h. The calculation transforms each replicate's \(k_e\) first; it does not transform the mean or standard error of \(k_e\) to obtain this result.

Why Joint Parameter Dependence Matters

A derived quantity often depends on more than one parameter. The earlier material on Bootstrap Covariance and Parameter Correlation describes how bootstrap estimates can covary. For a derived quantity, that dependence is not incidental: it can determine the shape and spread of the transformed distribution. Pairing means that all parameters used in \(g\) come from the same replicate.

Worked Example: Clearance and Exposure

In four hypothetical bootstrap fits, the estimated clearances are 8, 10, 12, and 10 L/h. For a fixed hypothetical IV dose of 100 mg under linear PK, exposure is calculated as \(\mathrm{AUC}=\mathrm{Dose}/CL\). The replicate AUC values are therefore \(100/8=12.50\), \(100/10=10.00\), \(100/12=8.333\), and \(100/10=10.00\) mg·h/L.

The mean of these four transformed values is \(40.833/4=10.208\) mg·h/L. Their sample standard deviation is:

s_{\mathrm{AUC}} = \sqrt{\frac{ (12.50-10.208)^2+(10.00-10.208)^2+ (8.333-10.208)^2+(10.00-10.208)^2 }{3}} \approx 1.718\ \mathrm{mg\cdot h/L}.

This standard deviation describes the spread of exposure calculated from these clearance replicates. In an actual analysis, many more successful replicates would be used, and an interval would be based on the transformed replicate distribution. The dose is fixed here; if dose were also estimated or uncertain, its replicate value would need to remain paired with clearance when calculating AUC.

The same principle applies when the quantity depends on correlated parameters. For instance, in a one-compartment model with first-order elimination, \(t_{1/2}=\ln(2)V/CL\). The clearance and volume in that expression should be taken from the same bootstrap fit. Recombining them across fits destroys the pairing and can change the derived distribution even when the separate marginal distributions of clearance and volume remain unchanged.

Worked Example: Half-Life From Paired Clearance and Volume Estimates

Consider four hypothetical bootstrap fits from a one-compartment model. The estimated pairs \((CL,V)\), with \(CL\) in L/h and \(V\) in L, are \((8,40)\), \((10,50)\), \((12,72)\), and \((10,60)\). Assume the elimination half-life formula applies to this model.

The half-lives calculated within the four replicates are \(0.6931(40/8)=3.466\), \(0.6931(50/10)=3.466\), \(0.6931(72/12)=4.159\), and \(0.6931(60/10)=4.159\) h. Their mean is approximately \(3.812\) h, and their sample standard deviation is:

s_{t_{1/2}} = \sqrt{\frac{2(3.466-3.812)^2+2(4.159-3.812)^2}{3}} \approx 0.400\ \mathrm{h}.

Now keep the same four clearance estimates and the same four volume estimates, but incorrectly pair the volumes in the order 72, 60, 40, 50 L. The resulting half-lives would be approximately 6.238, 4.159, 2.310, and 3.466 h. The marginal lists of clearance and volume have not changed, but the calculated half-life distribution has. This demonstrates why parameter estimates must be transformed as replicate-specific vectors, not as separate lists or marginal summaries.

Justification and the Delta-Method Comparison

The bootstrap transformation is direct: if the replicate vectors represent the sampling variability of the fitted parameter vector, applying \(g\) to each vector carries that variability through the function. A local approximation explains why covariance matters and when a simpler alternative may be adequate. For a differentiable \(g\), a first-order Taylor expansion around an estimate \(\hat{\boldsymbol{\theta}}\) gives:

g(\hat{\boldsymbol{\theta}}+\boldsymbol{\delta}) \approx g(\hat{\boldsymbol{\theta}}) + \nabla g(\hat{\boldsymbol{\theta}})^{T}\boldsymbol{\delta}.

If the parameter estimate has covariance matrix \(\Sigma\), the linear term has approximate variance \(\nabla g^{T}\Sigma\nabla g\). For \(g(V,CL)=\ln(2)V/CL\), this approximation uses both variances and the covariance of \(V\) and \(CL\); omitting the covariance term is generally unjustified. The bootstrap instead evaluates the nonlinear function across the observed joint replicate distribution. It can therefore reflect nonlinear effects and asymmetry that a first-order approximation may miss, provided the bootstrap distribution itself is informative.

The delta method can be useful when the function is smooth near the estimate and parameter uncertainty is small enough for a local linear approximation. It can be less reliable for strongly nonlinear transformations, ratios with uncertain denominators, or skewed distributions. Direct bootstrap propagation is not automatically reliable in those cases either: it depends on adequate, representative bootstrap replicates and on the appropriateness of the fitted model and resampling procedure.

Diagnostics, Failure Modes, and Limitations

Inspect the transformed replicates, not only their standard error or interval. Check the range, histogram or quantiles, and whether values are finite and physiologically interpretable. For a positive quantity such as clearance or half-life, negative transformed values can indicate a formula or parameterization error. Extremely large ratios may indicate a denominator close to zero; such replicates should not be silently discarded merely because they are inconvenient.

Check the formula and units. For example, dose in mg divided by clearance in L/h gives AUC in mg·h/L. A unit mismatch can make an otherwise correct transformation meaningless.
Keep parameter vectors intact. Sorting, independently resampling, or pairing estimates from different replicates can erase joint dependence.
Account for failed fits. The transformed distribution is based on successful fits, so a high or selective failure rate can make it unrepresentative. Report and investigate failures under the bootstrap workflow rather than treating success as automatic.
Check Monte Carlo stability. Repeat summaries or compare batches of replicates. A result that changes materially as more replicates are added may be too imprecise for the intended interpretation.
Check the model assumptions behind the quantity. The one-compartment half-life formula is not interchangeable with a terminal half-life formula for a multicompartment model. Use a formula that matches the model and the clinical definition of the target.

A bootstrap distribution quantifies uncertainty under its resampling scheme and fitted model; it does not correct model misspecification or create information absent from the data. Derived quantities may also be more weakly identified than their component parameters. A seemingly stable parameter estimate does not guarantee a stable ratio, especially when the denominator has substantial uncertainty.

Common pitfall: “I will calculate the half-life from the mean volume and mean clearance, then combine their standard errors.” That generally does not give the bootstrap distribution or uncertainty of half-life. The correction is to calculate the formula from each replicate's paired volume and clearance estimates, then summarize those transformed values.

Check Your Understanding

Answer each question, then compare your reasoning with the short solutions below.

Each bootstrap fit provides a clearance estimate and a volume estimate. What must be done before calculating a bootstrap standard error for \(t_{1/2}=\ln(2)V/CL\)?
A hypothetical replicate has dose 150 mg and clearance 15 L/h. Calculate AUC and give its units.
Why can independently pairing clearance and volume estimates change the bootstrap uncertainty of half-life, even if the two marginal lists are unchanged?
Name one situation in which the delta method may be less reliable than direct bootstrap propagation.

Solutions. (1) Calculate \(\ln(2)V/CL\) within every replicate using that replicate's paired estimates, then summarize the transformed values. (2) \(\mathrm{AUC}=150/15=10.00\) mg·h/L. (3) The derived quantity depends on both parameters, so their within-replicate dependence affects the ratios and their distribution; changing pairings changes those ratios. (4) One example is a strongly nonlinear ratio with an uncertain denominator, for which a local linear approximation may not represent the transformed distribution well.