Tutorials › Nonlinear Mixed-Effects Models (NLME) › A Worked Subject Bootstrap for Clearance Uncertainty

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

A Worked Subject Bootstrap for Clearance Uncertainty

Follow a subject-level bootstrap from longitudinal concentration data through model refitting to percentile and BCa uncertainty intervals for typical clearance.

Advanced 12 min read

What You'll Learn

  • Resample subjects as clusters so each subject’s full concentration profile remains intact.
  • Define the target as typical clearance in a one-compartment NLME model.
  • Calculate percentile interval endpoints from ordered subject-bootstrap estimates.
  • Calculate BCa bias correction and acceleration, then adjust the interval probabilities.
  • Diagnose failed refits, influential subjects, and unstable bootstrap tails.
  • Compare subject bootstrap intervals with Wald, profile-likelihood, and parametric-bootstrap approaches.

Why Resample Subjects?

A confidence interval for typical clearance should reflect uncertainty in the population of subjects, not pretend that every concentration measurement came from an independent person. In a longitudinal study, measurements within a subject share that subject’s latent pharmacokinetic characteristics and are generally dependent. A subject bootstrap preserves this structure: it resamples whole subjects and carries all of each selected subject’s records into the refitted data set.

The worked example uses a small, hypothetical intravenous pharmacokinetic study. Its purpose is to make the resampling and interval calculations explicit. With only six subjects, the resulting intervals should be viewed as an illustration of the method, not as a claim that a small sample guarantees reliable tail estimates.

Definition: In a subject bootstrap, each bootstrap replicate is formed by sampling subjects with replacement from the observed subjects. All records belonging to a sampled subject are included together. The model is refitted to each replicate, producing a bootstrap estimate of the parameter of interest.

Model and Resampling Assumptions

Suppose six subjects each receive a 100 mg intravenous bolus, and concentration is measured at 0.5, 2, and 6 hours. Consider a one-compartment model with a shared volume of distribution \(V\) and subject-specific clearance:

$$ C_{ij}=\frac{100\ \mathrm{mg}}{V}\exp\left(-\frac{CL_i}{V}t_{ij}\right)+\epsilon_{ij}, \qquad CL_i=CL_{\mathrm{typ}}\exp(\eta_i), \qquad \eta_i\sim N(0,\omega_{CL}^{2}). $$

Here \(C_{ij}\) is the concentration for subject \(i\) at time \(t_{ij}\), \(CL_{\mathrm{typ}}\) is typical clearance in L/h, and \(\eta_i\) is the subject-level random effect. For this illustration, assume independent residual errors with a common normal distribution, independent subject random effects, and independent subjects sampled from the target population. In an actual analysis, the residual-error model and random-effects structure must be chosen and checked for the data; the bootstrap does not make a misspecified model correct.

The subject bootstrap treats the observed subjects as the empirical sample from the target population. It assumes that subjects are the independent sampling units and that the study’s sampling scheme is reasonably represented by those subjects. It is not appropriate to resample individual concentration records independently: that would break within-subject profiles and understate the role of between-subject variability. If subjects were sampled in clusters, such as multiple subjects per clinic, the independent resampling unit may instead need to be the clinic.

Worked Example: Assemble the Longitudinal Data

The following hypothetical concentrations are in mg/L. Each row is one subject, and the three measurements in that row are the complete observed profile used in the fit.

Subject0.5 h2 h6 h
A4.473.081.13
B4.613.341.47
C4.322.760.89
D4.733.421.54
E4.503.141.20
F4.392.931.01

For the illustrative analysis, suppose the model is fitted to all six profiles and gives \(\widehat{CL}_{\mathrm{typ}}=6.20\) L/h. This is a hypothetical fit result, supplied to demonstrate uncertainty calculations; it is not obtained by taking an arithmetic average of the table. The model estimates its parameters by using the full concentration-time profiles and the specified random-effects and residual-error distributions.

A bootstrap replicate might select subjects A, A, C, D, F, F. It contains six subject slots, but only four distinct original profiles. Both copies of A retain all three A measurements and both copies of F retain all three F measurements. The copies are treated as separate sampled subject instances in the refit; their records are not split across the data set.

From Replicates to a Percentile Interval

For each replicate, refit the same model with the same parameterization, residual-error assumptions, estimation method, and constraints. Record the resulting \(\widehat{CL}_{\mathrm{typ}}^{*(b)}\), where \(b\) indexes the bootstrap replicate. The empirical distribution of these estimates approximates the sampling distribution under repeated sampling of subjects, subject to the model and bootstrap assumptions.

For \(B\) estimates sorted from smallest to largest, the two-sided 95% percentile interval uses the empirical 2.5th and 97.5th percentiles. Its endpoints are quantiles of the bootstrap estimates, not the original estimate plus or minus a standard error. The precise interpolation or rank convention should be documented, particularly when \(B\) is modest.

Worked Example: Percentile Endpoints

Suppose \(B=2{,}000\) subject-bootstrap refits all converge. In this hypothetical output, the estimates at the rounded ranks nearest \(2{,}000(0.025)=50\) and \(2{,}000(0.975)=1{,}950\) are 3.62 and 8.76 L/h, respectively.

The percentile interval is therefore \([3.62,\ 8.76]\) L/h. The lower endpoint is the 2.5th percentile of the replicate estimates and the upper endpoint is their 97.5th percentile, according to the stated rank convention. In repeated use of this procedure under its assumptions, the interval method is intended to cover the fixed typical clearance at approximately 95%; this particular interval does not mean there is a 95% probability that the fixed parameter lies between its endpoints.

Bias-Corrected and Accelerated (BCa) Interval

A percentile interval uses the nominal tail probabilities directly. A BCa interval modifies them to account for estimated bias and acceleration, a measure of how quickly the estimator’s sampling variability changes across the parameter space. The adjustment uses two quantities calculated from the bootstrap and delete-one-subject fits.

Definition: Let \(\widehat{\theta}\) be the estimate from the original data and \(\widehat{\theta}^{*(b)}\) the estimate from bootstrap replicate \(b\). The BCa bias-correction term is \(z_0=\Phi^{-1}\{P^*(\widehat{\theta}^{*}<\widehat{\theta})\}\), estimated by the fraction of bootstrap estimates below the original estimate. If \(\widehat{\theta}_{(i)}\) is the estimate after deleting subject \(i\), define \(\bar{\theta}_{(.)}\) as the mean of the delete-one estimates. The acceleration is $$ a=\frac{\sum_{i=1}^{n}\left(\bar{\theta}_{(.)}-\widehat{\theta}_{(i)}\right)^3} {6\left[\sum_{i=1}^{n}\left(\bar{\theta}_{(.)}-\widehat{\theta}_{(i)}\right)^2\right]^{3/2}}. $$

For a nominal tail probability \(\alpha\), its BCa-adjusted probability is

$$ \alpha_{\mathrm{BCa}}= \Phi\left[ z_0+ \frac{z_0+z_{\alpha}} {1-a(z_0+z_{\alpha})} \right], \qquad z_{\alpha}=\Phi^{-1}(\alpha). $$

The lower and upper interval probabilities are calculated separately, with \(\alpha=0.025\) and \(\alpha=0.975\). The corresponding quantiles of the bootstrap estimates give the BCa endpoints. When the denominator is near zero, the adjustment can become extreme; that is a warning about unstable acceleration or a problematic bootstrap distribution, not a reason to force an endpoint.

Calculating the BCa Adjustment

The bias correction shifts the nominal quantiles according to where the original estimate falls in the bootstrap distribution. Acceleration adjusts for asymmetry in the estimator’s sampling behavior as different subjects are omitted. The following calculations use hypothetical summary output for the same six-subject example.

Worked Example: Bias Correction and Acceleration

Suppose 1,040 of the 2,000 bootstrap estimates are below the original fit, 6.20 L/h. Then

$$ z_0=\Phi^{-1}\left(\frac{1{,}040}{2{,}000}\right) =\Phi^{-1}(0.5200)\approx 0.0502. $$

Suppose the six delete-one fits give \(\sum_i(\bar{\theta}_{(.)}-\widehat{\theta}_{(i)})^2=1.54\) and \(\sum_i(\bar{\theta}_{(.)}-\widehat{\theta}_{(i)})^3=0.168\), in the corresponding powers of L/h. Then

$$ a=\frac{0.168}{6(1.54)^{3/2}} \approx\frac{0.168}{11.47} \approx 0.0147. $$

For the lower endpoint, \(z_{0.025}\approx-1.960\). Substitution gives

$$ \Phi\left(0.0502+ \frac{0.0502-1.960}{1-0.0147(0.0502-1.960)}\right) =\Phi(-1.808)\approx 0.0353. $$

For the upper endpoint, \(z_{0.975}\approx1.960\), so

$$ \Phi\left(0.0502+ \frac{0.0502+1.960}{1-0.0147(0.0502+1.960)}\right) \approx\Phi(2.122)\approx0.9831. $$

Using the nearest-integer rank convention, the adjusted ranks are \(2{,}000(0.0353)=70.6\), rounded to 71, and \(2{,}000(0.9831)=1{,}966.2\), rounded to 1,966. Suppose the sorted bootstrap estimates at those ranks are 3.48 and 9.08 L/h. The BCa interval is \([3.48,\ 9.08]\) L/h. Its endpoints are wider than the percentile interval in this hypothetical output because the adjustment selects more extreme empirical quantiles.

Justification and Practical Checks

The subject-level resampling logic can be stated directly. If the observed independent subjects are the sampling units, resampling \(n\) subjects with replacement approximates drawing a new sample of size \(n\) from the empirical distribution of observed subject profiles. Refitting in each draw lets the estimate respond to changes in both the collection of profiles and the nonlinear estimation problem. Resampling measurements separately would instead approximate a different sampling process in which records, rather than subjects, were independent.

BCa’s adjustment follows from correcting the nominal normal-scale tail \(z_\alpha\) first for bias and then for acceleration. With \(a=0\), the denominator in the adjustment is 1, so only the bias shift \(z_0\) remains. With \(z_0=0\) as well, \(\alpha_{\mathrm{BCa}}=\Phi(z_\alpha)=\alpha\), and the BCa quantiles reduce to the percentile quantiles. This limiting case checks the direction and structure of the formula.

Before interpreting either interval, inspect the refit process and the replicate distribution. Check that the intended number of replicates converged, that estimates are finite and biologically plausible, and that boundary solutions or extreme values are understood. Examine the distribution for skewness, multimodality, long tails, or clusters associated with repeated copies of a small number of subjects. Compare the original estimate with the bootstrap distribution and inspect leave-one-subject-out estimates, since one influential profile can strongly affect both the estimate and the BCa acceleration.

If refits fail, report how many failed and investigate the reason. Repeated convergence failures may reveal weak identification, poor initial values, or a model that is unstable for plausible subject samples. Silently discarding failures can change the empirical distribution and make the interval misleading. Increasing \(B\) reduces Monte Carlo variability in estimated quantiles, but it does not repair model misspecification, a nonrepresentative sample, or a very small number of independent subjects.

How This Approach Compares and Where It Is Limited

ApproachWhat it usesImportant qualification
Wald intervalEstimate and an estimated standard error, often with a normal approximationCan be convenient, but may be inaccurate when the sampling distribution is skewed or near a parameter boundary.
Subject-bootstrap percentile or BCa intervalRepeated refits to resampled whole-subject dataReflects empirical subject variation but depends on adequate independent subjects and stable refits.
Profile-likelihood intervalLikelihood changes as the target parameter is fixed and nuisance parameters are reoptimizedUses likelihood curvature and optimization rather than empirical resampling; it can be informative when a Wald approximation is poor.
Parametric bootstrapNew data simulated from the fitted model, followed by refittingCan generate more varied profiles than resampling a tiny observed sample, but relies strongly on the assumed data-generating model.

These methods answer the same uncertainty question through different approximations; agreement is reassuring but not guaranteed. A subject bootstrap cannot create information absent from the observed subjects. In particular, a small sample may yield a discrete or unstable empirical tail, and a few distinct subjects can dominate many replicates. BCa is not automatically superior: acceleration based on delete-one estimates can be unstable when deleting one subject changes identifiability or causes a different convergence outcome. Report the target parameter, model, resampling unit, number of replicates, endpoint convention, convergence handling, and interval method so the calculation can be understood and reproduced.

Key takeaway: For longitudinal NLME data, resample independent subjects with their complete profiles, refit the same model in every replicate, and inspect the resulting estimates before interpreting percentile or BCa limits. The interval describes uncertainty under the stated sampling and model assumptions, not uncertainty about whether the model is correct.

Check Your Understanding

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

  1. Why should the three concentration records from one subject be resampled together rather than separately?
  2. In the hypothetical bootstrap output, what does the percentile interval \([3.62,\ 8.76]\) L/h mean?
  3. If \(z_0=0\) and \(a=0\), what happens to the BCa adjusted probabilities?
  4. What should be investigated if many bootstrap refits fail to converge?
  5. Does increasing the number of bootstrap replicates correct bias caused by a misspecified residual-error model?

Solutions. (1) The records within a subject form a longitudinal profile and are not independent sampling units; keeping them together preserves their subject-level dependence. (2) It is the empirical 2.5th-to-97.5th percentile interval of the replicate estimates under the chosen subject-bootstrap procedure, with endpoints in L/h; it is not a 95% probability statement about the fixed parameter. (3) The formula gives \(\alpha_{\mathrm{BCa}}=\Phi(z_\alpha)=\alpha\), so the adjusted probabilities equal the nominal ones. (4) Check model identifiability, initialization, boundaries, data composition, and whether failures are associated with particular resampled subjects; report the failures rather than silently omitting them. (5) No. More replicates reduce Monte Carlo error in the bootstrap calculation but do not correct misspecification or other systematic problems.