Why Simulate and Refit?
An NLME parameter estimate depends on the particular subjects and observations in a study. A Wald interval summarizes uncertainty using local curvature of the fitted likelihood; that summary can be unreliable when the likelihood is asymmetric, a parameter is near a boundary, or the model is strongly nonlinear. A parametric bootstrap offers another way to assess uncertainty: use the fitted model to generate new datasets, then apply the same estimation procedure to each one.
The result is not a collection of new observations from the actual study. It is a collection of hypothetical datasets generated under the fitted model and the study’s design. Consequently, the bootstrap describes how the estimator varies if the assumed model and its fitted parameters are a reasonable approximation to the data-generating process.
Model Assumptions and the Bootstrap Algorithm
For clarity, suppose subject \(i\) has observations \(y_{ij}\) at the recorded design points, such as dose and sampling times. A general hierarchical representation is
Here \(f\) is the structural model, \(\theta\) contains population parameters, \(\eta_i\) is the subject-specific random effect, and \(\Omega\) is its covariance matrix. The residual errors \(\epsilon_{ij}\) follow the residual model specified for the analysis; for example, they may be independent normal errors with mean zero and variance \(\sigma^2\). Other residual structures require simulation from those structures, not automatic use of this simple normal form.
A model-based bootstrap relies on the fitted model being an adequate description of the relevant data-generating features: the structural relationship, random-effects distribution, residual distribution and dependence, and any covariate relationships represented in the model. It also requires the fitted model to be estimable well enough that refitting simulated datasets is meaningful. A bootstrap does not remove these assumptions; it uses them to generate data.
Fit the specified NLME model to the original data and record the estimates \(\hat{\theta}\), \(\hat{\Omega}\), and residual parameters, together with the study design.
For each subject, simulate random effects from the fitted random-effects distribution, then simulate observations under the fitted model and residual distribution at the recorded design points.
Apply the same model, estimation method, parameter constraints, and relevant fitting settings to the replicate. Save the parameter estimates and record fitting diagnostics.
Repeat for \(B\) replicates. Use the resulting estimates to examine their sampling variation and, with a stated interval method, quantify uncertainty.
The essential idea is that the estimator is applied again, not merely that outcomes are simulated. If \(\hat{\psi}\) denotes the original estimate of a parameter vector \(\psi\), and \(\hat{\psi}^{*(b)}\) is the estimate from replicate \(b\), the bootstrap distribution is the collection \(\hat{\psi}^{*(1)},\ldots,\hat{\psi}^{*(B)}\). It approximates the repeated-sampling behavior of the estimator under the fitted model.
Worked Examples
Worked Example: Simulating Concentrations at Recorded Times
Hypothetical setup. Consider a one-compartment IV bolus model with dose \(D=100\) mg, volume \(V=20\) L, and individual clearance \(CL_i=4.421\exp(\eta_i)\) L/h, where \(\eta_i\sim N(0,0.09)\). The residual error is additive normal with standard deviation \(0.10\) mg/L. The fitted typical clearance is \(4.421\) L/h. For one illustrative simulated subject, take a random-effect draw displayed as \(\eta_i=0.000\), and residual draws \(+0.10\) mg/L at 1 hour and \(-0.05\) mg/L at 4 hours. These are hypothetical values, not study results.
The predicted concentration for this subject is \(D/V\exp(-CL_i t/V)\). At 1 hour, \(CL_i=4.421\exp(0.000)=4.421\) L/h, so
Adding the simulated residual gives \(4.008+0.10=4.108\) mg/L. At 4 hours,
so the simulated observation is \(2.0652-0.05=2.0152\) mg/L, or approximately \(2.015\) mg/L. In an actual replicate, this simulation is done for all subjects and their recorded design points; the replicate dataset is then refitted.
Worked Example: Calculating a Bootstrap Standard Error
Suppose six successful hypothetical refits give estimates \(0.01,\ 0.04,\ 0.08,\ 0.12,\ 0.15,\ 0.25\) for a log-clearance parameter. These deliberately small numbers illustrate the arithmetic; they are not enough replicates for a reliable final analysis. The mean of these estimates is
The sum of their squares is \(0.1075\). The sum of squared deviations from their mean can be calculated as
The sample standard deviation of the replicate estimates, used here as the bootstrap estimate of standard error, is therefore
Its units are squared-log-clearance units under the variance calculation, and log-clearance units for the standard error. It estimates the variability of the estimator under the fitted model; it is not the residual standard deviation of concentrations. In practice, substantially more replicates are needed to characterize uncertainty well.
Worked Example: Reading an Illustrative Bootstrap Interval
Suppose a fitted model estimates typical clearance as \(4.421\) L/h, so its log-scale estimate is \(\hat{\beta}=\log(4.421)\approx1.4864\). In a hypothetical run of \(B=2{,}000\) successful refits, suppose the 2.5th and 97.5th percentiles of the simulated log-clearance estimates are 1.31 and 1.67. A central empirical bootstrap interval on the log scale is thus \([1.31,1.67]\). Since exponentiation is increasing, transforming both endpoints gives
This is an illustrative 95% percentile bootstrap interval: its endpoints are taken from the tails of the simulated estimate distribution, then expressed on the clearance scale. It is not a probability statement that the fixed clearance has a 95% chance of lying in this particular interval. Its repeated-sampling interpretation depends on the model and simulation scheme being adequate. The detailed properties and alternatives for bootstrap interval construction require separate consideration.
Why the Plug-In Procedure Is Reasonable—and Where It Can Fail
The justification is an approximation to repeated sampling. If the fitted model were the data-generating model, repeatedly collecting studies under that model and refitting would produce a sampling distribution for the estimator. The true parameter values are unknown, so the parametric bootstrap substitutes the original estimates for them and simulates from that fitted model. Repeating the fitting step then approximates the distribution of the estimator under this plug-in data-generating process.
This approximation is most defensible when the model is well specified for the feature of interest and the estimator behaves regularly in the region of the fitted parameter values. It may be poor if the structural model misses systematic features, the random-effects distribution is seriously wrong, important residual correlation is omitted, or parameters are weakly identified. In those cases the simulated datasets reproduce the fitted assumptions, not the data features the model failed to represent.
Inspect the replicate estimates and fitting records rather than relying only on a reported interval. Useful checks include whether estimates are systematically shifted relative to the original estimate, whether distributions are strongly skewed or have multiple clusters, whether variance estimates accumulate near zero, and whether simulated datasets commonly produce implausible predictions. Compare observed data features with simulated ones using appropriate predictive checks. Frequent convergence warnings or boundary estimates are also evidence that the fitted procedure may be unstable; they should not be hidden by summarizing only convenient output.
A failed fit matters because the bootstrap distribution is meant to reflect the estimator actually used. If some replicates cannot be fit, simply discarding them can change the distribution being summarized. The reasons for failures and their consequences need explicit investigation; a numerical interval based on a selectively retained subset can be misleading.
How It Compares With Other Approaches
| Approach | What generates uncertainty information | Main consideration |
|---|---|---|
| Parametric bootstrap | Simulated datasets from the fitted NLME model, each refit using the analysis procedure | Can reflect nonlinear refitting behavior, but depends on the fitted data-generating assumptions |
| Subject-level nonparametric bootstrap | Resampled subjects, with their observed records, followed by refitting | Relies on the observed subjects representing the population; resampling cannot generate features absent from the sample |
| Wald interval | Local curvature of the fitted likelihood, summarized by an estimated standard error | Computationally direct, but local symmetry and approximation conditions may be inadequate |
These approaches answer related but not identical questions. A parametric bootstrap is attractive when the model gives a credible account of the variability and the observed sample is too limited to represent every relevant pattern directly. It is not automatically superior to a subject-level bootstrap or a likelihood-based interval: model adequacy, estimability, and the parameter’s sampling behavior should guide the choice.
Common Pitfall and Practical Limits
Another mistake is to simulate only residual errors around fixed individual predictions and then call the result a population-level parametric bootstrap. That procedure omits between-subject variability if the target simulation requires new random effects. Conversely, which random effects are redrawn and which quantities are held fixed defines the simulation target; the conditional and unconditional choices are treated separately in the relevant material. State the target clearly rather than treating every simulation scheme as interchangeable.
Finally, Monte Carlo variability remains even with a correct model: a finite number \(B\) of replicates gives an imperfect approximation to the bootstrap distribution. A large \(B\) cannot repair model misspecification or poor identifiability, while a small \(B\) can make tail estimates unstable. The parametric bootstrap is best viewed as a model-dependent uncertainty tool, accompanied by diagnostics and transparent reporting of the simulation and refitting procedure.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- In an NLME parametric bootstrap, what varies across replicate datasets, and what is refitted?
- Why does simulating outcomes from the fitted model not demonstrate that the model fits the observed data?
- For bootstrap estimates \(0.01,0.04,0.08,0.12,0.15,0.25\), what is their mean and sample standard deviation?
- If log-scale interval endpoints are 1.31 and 1.67 for log clearance, what interval do they imply on the clearance scale?
Solutions. (1) Simulated subject effects and residual errors vary according to the stated fitted model; the same NLME analysis is refitted to each generated dataset. (2) The simulation assumes the fitted model and therefore cannot independently validate it; observed and simulated features must be compared. (3) The mean is \(0.65/6=0.1083\). The squared-deviation sum is \(0.1075-0.65^2/6=0.0370833\), so the sample standard deviation is \(\sqrt{0.0370833/5}=0.08612\). (4) Exponentiating the endpoints gives \([\exp(1.31),\exp(1.67)]\approx[3.706,5.312]\) L/h.