Tutorials › Nonlinear Mixed-Effects Models (NLME) › Percentile Bootstrap Confidence Intervals in NLME Models

Bootstrap, Profile Likelihood and Parameter Uncertainty · Tutorial 622 of 797

Percentile Bootstrap Confidence Intervals in NLME Models

You will calculate percentile intervals for model and derived parameters and learn why a skewed bootstrap distribution alone does not guarantee nominal coverage.

Advanced 10 min read

What You'll Learn

  • Define percentile bootstrap confidence limits using empirical quantiles of refitted parameter estimates.
  • Calculate interpolated endpoints using a stated quantile convention.
  • Form intervals for positive and derived NLME parameters.
  • Explain the coverage condition and why skewness can invalidate nominal coverage.
  • Identify diagnostics and limitations that affect interpretation of an interval.

What a Percentile Interval Uses

A bootstrap analysis produces a collection of parameter estimates by repeatedly resampling or simulating data under a specified design and refitting the NLME model. The percentile method turns that collection into an interval by taking its lower and upper empirical quantiles. This is attractive when the estimates are skewed: unlike a symmetric Wald interval, the two distances from the fitted estimate to the endpoints need not be equal.

The method’s calculation is simple, but its coverage is not automatic. In particular, reproducing a skewed sampling distribution in bootstrap replicates does not by itself establish that the resulting interval covers the true parameter at its nominal rate. The distinction between a distributional approximation and a coverage argument is important when interpreting NLME parameter uncertainty.

Definition: Let \(\hat{\theta}\) be the estimate from the original data, and let \(\hat{\theta}^{*(1)},\ldots,\hat{\theta}^{*(B)}\) be the estimates from \(B\) valid bootstrap replicates. The percentile \(100(1-\alpha)\%\) interval for \(\theta\) is the interval from the \(\alpha/2\) to the \(1-\alpha/2\) quantile of the bootstrap estimates. Its endpoints are empirical quantiles of \(\hat{\theta}^{*}\), not necessarily \(\hat{\theta}\) plus or minus the same distance.

The parameter \(\theta\) may be a structural parameter, a variance component, or a derived quantity computed from each replicate. The interval’s units are the units of that parameter. For a 95% interval, the target quantile probabilities are 0.025 and 0.975; other confidence levels use the corresponding tail probabilities.

Before interpreting the endpoints, the bootstrap must represent the relevant sampling process. Resampling must use the appropriate unit, preserve any relevant study strata, and retain the dosing and sampling design. A parametric bootstrap must simulate from a model and parameter values that are defensible for the purpose. Replicate fits must also yield estimates of the same parameter on a comparable scale. These conditions do not guarantee good coverage, but violations undermine the bootstrap distribution before the percentile calculation even begins.

Calculating Empirical Quantiles

An empirical quantile depends on the convention used to locate it between ordered estimates. Software packages can use different conventions, so a reproducible analysis should identify the convention or software used. In the examples below, use linear interpolation with position \(h=1+(B-1)p\) for quantile probability \(p\). If \(h=j+\gamma\), where \(j\) is an integer and \(0\leq\gamma\leq1\), interpolate as \((1-\gamma)\hat{\theta}^{*}_{(j)}+\gamma\hat{\theta}^{*}_{(j+1)}\), where parentheses denote order from smallest to largest.

Worked Example: A Clearance Interval With a Longer Upper Tail

Suppose 11 hypothetical valid bootstrap fits give the following ordered estimates of clearance divided by bioavailability, in L/h: 3.8, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.7, 5.0, 5.4, and 6.2. Calculate a 90% percentile interval.

For a 90% interval, \(\alpha=0.10\), so the endpoint probabilities are 0.05 and 0.95. Here \(B=11\), and the interpolation position is \(h=1+(11-1)p\). At \(p=0.05\), \(h=1.5\): the lower endpoint is halfway between the first two estimates, \(3.8+0.5(4.0-3.8)=3.9\) L/h. At \(p=0.95\), \(h=10.5\): the upper endpoint is halfway between the tenth and eleventh estimates, \(5.4+0.5(6.2-5.4)=5.8\) L/h.

The resulting 90% percentile interval is 3.9 to 5.8 L/h. Its upper side is farther from the original-data estimate if that estimate is near the centre of the bootstrap distribution, reflecting the long upper tail of these illustrative replicates. This describes the calculated interval; it does not, by itself, demonstrate that repeated use of the procedure would cover the true clearance with probability 0.90.

With a finite number of replicates, empirical quantiles have Monte Carlo variability. A small \(B\) can make tail quantiles especially unstable, and interpolation cannot create information beyond the observed replicates. The interval’s reported precision should therefore be sensible relative to the number of bootstrap fits and the variability of its endpoints.

Derived Parameters and Scale

For a derived pharmacokinetic quantity, calculate that quantity separately in every bootstrap replicate, then take empirical quantiles of the derived estimates. This preserves the replicate-level relationship between the model parameters that determine the derived quantity. For a one-to-one monotone transformation, transforming endpoints can sometimes produce the same set of quantiles with the order reversed when the transformation is decreasing. But the direction must be handled correctly, and the direct replicate-by-replicate calculation is a reliable way to avoid mistakes.

Worked Example: A Half-Life Interval From Bootstrap Rates

Suppose 10 hypothetical bootstrap fits produce elimination-rate estimates, ordered from smallest to largest, of 0.055, 0.060, 0.064, 0.068, 0.072, 0.078, 0.085, 0.095, 0.110, and 0.140 h\(^{-1}\). For each fit, calculate \(t_{1/2}=\ln(2)/k\), in hours. The resulting half-lives, now ordered from smallest to largest, are approximately 4.951, 6.301, 7.296, 8.155, 8.886, 9.627, 10.193, 10.830, 11.552, and 12.603 h. For instance, the largest rate gives \(0.693147/0.140=4.951\) h, while the smallest gives \(0.693147/0.055=12.603\) h.

Calculate an 80% percentile interval for half-life. The endpoint probabilities are 0.10 and 0.90. With \(B=10\), the positions are \(h=1+9(0.10)=1.9\) and \(h=1+9(0.90)=9.1\). Thus, the lower limit is \(4.951+0.9(6.301-4.951)=6.166\) h. The upper limit is \(11.552+0.1(12.603-11.552)=11.657\) h.

The illustrative 80% interval is 6.166 to 11.657 h. It is calculated from the derived half-lives, not by taking the 10th and 90th percentiles of the rate estimates in the same order: the transformation \(\ln(2)/k\) is decreasing, so small rates correspond to large half-lives. Ten replicates are far too few for stable tail estimation in practice; they are used here only to show the arithmetic.

Why Percentile Intervals May Miss Their Nominal Coverage

A percentile interval is often motivated by approximating the estimator’s sampling distribution with the bootstrap distribution. The coverage implication requires a further step. Write the estimation error as \(E=\hat{\theta}-\theta\), and let \(g_p\) be its \(p\) quantile. If the bootstrap distribution of \(\hat{\theta}^{*}-\hat{\theta}\) adequately approximates the sampling distribution of \(E\), then the percentile endpoints are approximately \(\hat{\theta}+g_{\alpha/2}\) and \(\hat{\theta}+g_{1-\alpha/2}\). For the true value to fall between them, the estimation error must satisfy the following event:

$$ -g_{1-\alpha/2}\ \leq\ E\ \leq\ -g_{\alpha/2}. $$

Approximating the bootstrap error distribution with the sampling error distribution alone does not make this event have probability \(1-\alpha\). A sufficient condition is that the error distribution is symmetric about zero, so \(g_p=-g_{1-p}\); under that condition, the event becomes \(g_{\alpha/2}\leq E\leq g_{1-\alpha/2}\), which has probability \(1-\alpha\) for a continuous error distribution. More generally, nominal coverage requires the probability of the displayed event to be \(1-\alpha\), which need not hold for a biased or asymmetric error distribution.

Worked Example: A Skewed Error Distribution Does Not Guarantee Coverage

Consider a hypothetical estimation error \(E=Z-1\), where \(Z\) has an exponential distribution with mean 1. This error has mean zero but is strongly right-skewed. Its quantile is \(g_p=-\ln(1-p)-1\). For a nominal 90% percentile interval, \(g_{0.05}=-\ln(0.95)-1\approx-0.9487\), and \(g_{0.95}=-\ln(0.05)-1\approx1.9957\).

Under an ideal bootstrap approximation to this same error distribution, coverage is the probability that \(-g_{0.95}\leq E\leq-g_{0.05}\), or approximately \(-1.9957\leq E\leq0.9487\). Since \(E=Z-1\geq-1\), the lower bound excludes no possible values, and the upper bound is equivalent to \(Z\leq1.9487\). Therefore the coverage is \(P(Z\leq1.9487)=1-\exp(-1.9487)\approx0.8575\), not 0.90.

This constructed example shows that even an exact match between the bootstrap and sampling error distributions need not yield nominal percentile coverage when the error is asymmetric. In a real NLME analysis, the error distribution is not known, so the example is a justification for caution rather than a correction factor to apply to an interval.

Diagnostics, Alternatives, and Limitations

Inspect the bootstrap estimates as a distribution, not just as two numbers. A histogram or quantile summary can reveal skewness, extreme estimates, piles at a parameter boundary, or several distinct modes. Check whether these patterns are scientifically plausible and whether they reflect genuine sampling uncertainty, numerical instability, or a weakly identified parameter. Also check that the same parameterization and units were used in every fit. A percentile calculation cannot repair an inappropriate resampling scheme, an unstable estimation procedure, or a model that does not identify the target parameter well.

Only valid fits should contribute estimates, but excluding failed fits is not automatically harmless: failures may be related to particular regions of the resampled data or parameter space. The previous material on failed bootstrap fits emphasizes why their pattern and handling matter. Report the number of attempted and valid replicates, assess whether failures are concentrated in a consequential way, and do not silently treat the surviving fits as an unselected sample from the intended bootstrap distribution.

Compared with a Wald interval, the percentile interval uses empirical endpoints and can be asymmetric; a Wald interval instead relies on an estimated standard error and a normal-based construction. This flexibility is useful when a normal approximation is visibly questionable, but it is not a general guarantee of better coverage. The percentile method can inherit bias and skewness in an unhelpful way, and a finite bootstrap can give noisy tail limits. Other interval methods make different adjustments or use likelihood information, but choosing among methods requires considering the parameter scale, regularity, fit stability, and the intended coverage claim.

Key takeaway: A percentile interval is formed from the empirical lower and upper quantiles of bootstrap parameter estimates. It can reflect asymmetry in those estimates, but nominal coverage requires more than a good-looking bootstrap distribution: the bootstrap must represent sampling error, and the relevant error-quantile relation must support the desired coverage.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. For a 95% percentile interval, which quantiles of the bootstrap estimates are used?
  2. With the interpolation rule \(h=1+(B-1)p\), find the 25th percentile for ordered estimates 2, 4, 8, and 10.
  3. Why is a bootstrap distribution that reproduces an asymmetric estimation-error distribution not, by itself, enough to establish nominal percentile coverage?
  4. For \(t_{1/2}=\ln(2)/k\), do the smallest rates correspond to the smallest or largest half-lives?

Solutions. (1) Use the 0.025 and 0.975 quantiles. (2) Here \(B=4\), so \(h=1+3(0.25)=1.75\); the result is \(2+0.75(4-2)=3.5\). (3) Coverage depends on whether the estimation error falls between the negatives of the upper and lower error quantiles. For an asymmetric distribution, that event need not have probability \(1-\alpha\), even if the bootstrap reproduces the error distribution. (4) The smallest rates correspond to the largest half-lives because dividing \(\ln(2)\) by a smaller positive \(k\) gives a larger value.