Why the Search Itself Matters
A PK/PD analysis may compare several scientifically plausible exposure-response structures, covariate specifications, or residual models. Even when each candidate is fit correctly, selecting the best-looking result from many candidates changes what the reported fit and estimates mean. The selected model has won a competition: some of its advantage may reflect genuine structure, and some may be chance alignment between the candidate and the observed data.
This issue is not limited to formal hypothesis tests. It arises whenever the data influence which candidates are tried, which diagnostics trigger changes, which fit measure determines the winner, or which estimates are highlighted. A defensible candidate set, as discussed in the earlier material on constructing scientifically defensible PK/PD candidate sets, limits arbitrary search. It does not remove the uncertainty created when selection is based on the same data used for estimation.
Let \(M_1,\ldots,M_K\) denote a fixed set of candidate models, and let \(S_j\) be a fit score for model \(M_j\), with larger values defined to be better. The selected model is \(M_{\hat{j}}\), where \(\hat{j}=\arg\max_j S_j\). If the score can be decomposed as \(S_j=\mu_j+\varepsilon_j\), where \(\mu_j\) is the candidate's underlying expected score and \(\varepsilon_j\) is sampling variation, selection favors candidates with favorable realized errors as well as candidates with genuinely better expected scores.
Why Repeated Searches Bias Fit and Estimates
The basic mechanism can be seen when all candidates have the same underlying score, so \(\mu_j=\mu\), but their observed scores vary. Then the winning score is \(\mu+\max_j\varepsilon_j\). For integrable centered errors, \(E[\max_j\varepsilon_j]\ge 0\), because \(\max_j\varepsilon_j\ge\varepsilon_k\) for any \(k\) and \(E[\varepsilon_k]=0\); the inequality is strict unless all candidate errors are identical almost surely. Thus, the expected winning score is at least the underlying score, with equality when errors are perfectly shared; dependence can reduce or eliminate expected optimism.
Repeated testing gives a familiar special case. If \(m\) independent null hypotheses are each tested at significance level \(\alpha\), the probability of at least one false rejection is \(1-(1-\alpha)^m\). The independence assumption is essential for that exact formula. Without independence, the union bound gives a family-wise error probability no greater than \(m\alpha\), so a Bonferroni threshold of \(\alpha/m\) controls it at no more than \(\alpha\).
Worked Example: Chance of at Least One False Signal
Suppose a hypothetical analysis screens five exposure-response features, each by a separate test at \(\alpha=0.01\), and assume for calculation that the five tests are independent and all null hypotheses are true. What is the probability of at least one false rejection?
Each test has probability \(1-0.01=0.99\) of not rejecting falsely. Under independence, the probability that all five avoid a false rejection is \(0.99^5\). Thus,
The family-wise false-rejection probability is about 4.90%, not 1%. If independence is not credible, the exact calculation should not be used; the union bound instead gives an upper bound of \(5(0.01)=0.05\). A search over correlated models is not equivalent to five independent tests, but it still requires accounting for the fact that the winner was selected.
Parameter estimates can also be affected. Selection is based on data that also determine parameter estimates, so the reported estimate is conditional on having passed a data-dependent selection event. The direction and size of this effect depend on the selection rule, parameterization, and correlation between estimates and selection scores. It is not correct to assume that every parameter must be biased upward. A simple one-sided screening example makes the conditional effect explicit.
Suppose an investigator reports a positive treatment effect only if its standardized estimate \(Z=\hat{\theta}/SE\) exceeds 1.96. Under the null value \(\theta=0\), assume \(Z\sim N(0,1)\). Conditional on selection, the expected standardized estimate is the upper-tail truncated-normal mean:
where \(\phi\) and \(\Phi\) are the standard normal density and distribution functions. The integral identity follows because \(\phi'(z)=-z\phi(z)\), so \(\int_c^\infty z\phi(z)\,dz=\phi(c)\). At \(c=1.96\), the conditional mean is approximately \(0.05844/0.02500=2.338\), despite the true effect being zero. This example is a one-sided screen, not a universal formula for PK/PD selection. It shows why estimates from selected analyses can be systematically more extreme than estimates from a prespecified analysis.
Bootstrap the Entire Search Procedure
A bootstrap aimed at selection must reproduce the procedure that created the reported result, not merely refit the winning model. For repeated PK/PD observations, the observational unit for resampling is commonly the subject: sample subjects with replacement and retain each sampled subject's linked dose history, covariates, concentrations, and PD observations. This preserves within-subject dependence and the repeated-measures structure. Resampling individual records independently would generally break that structure and answer a different question.
Before resampling, specify the candidate set and freeze the search algorithm. This includes the data-processing rules, candidate fits, convergence and diagnostic rules, score or test used to compare models, stopping rules, and tie-breaking rule. If analysts would revise a model after a diagnostic, that revision rule belongs in the algorithm; silently omitting it from the bootstrap fails to represent the actual search. Resampling quantifies uncertainty for the specified procedure. It cannot make an unstable or scientifically inappropriate search valid.
State the target population, the subject as the resampling unit, the fixed candidate set, and the exact rule that selects a model.
For each bootstrap replicate, draw \(n\) subjects with replacement from the observed \(n\) subjects, carrying all records for each sampled subject into the replicate.
Refit every candidate and apply the same diagnostics, score, exclusions, convergence rules, and selection rule. Record the selected model and its parameter estimates.
Across replicates, report selection frequencies, parameter distributions with their conditioning made explicit, convergence failures, and—if estimating optimism—the corresponding training and evaluation scores.
If \(B\) bootstrap replicates are attempted and model \(M_j\) is selected in \(b_j\) successful replicates, its bootstrap selection frequency is \(b_j/B_{\text{success}}\), where \(B_{\text{success}}\) is the number of replicates in which the search completed under the prespecified rules. These frequencies describe the stability of selection under the empirical resampling scheme. They are not posterior probabilities that a model is true. A percentile interval for a parameter estimated only when \(M_j\) is selected is likewise a distribution conditional on that selection event; it is not automatically a confidence interval with guaranteed post-selection coverage.
Worked Example: Selection Frequencies From a Full-Search Bootstrap
In a hypothetical PK/PD analysis, 40 subjects contribute repeated concentration and response data. The prespecified search compares a linear exposure-response model, an Emax model, and an indirect-response model, selecting the model with the smallest AICc. A subject-level bootstrap runs \(B=200\) replicates and successfully completes all three fits in 200 replicates. The selected models are: Emax in 118 replicates, linear in 56, and indirect response in 26.
The estimated selection frequencies are
These frequencies sum to \((118+56+26)/200=1.000\). Although Emax is the most frequent winner, it is selected in only 59.0% of replicates. That is evidence of selection uncertainty: plausible resamples often favor another candidate. It does not mean Emax has a 59% probability of being biologically correct.
Suppose the full-data Emax fit estimates \(EC_{50}=1.20\) mg/L. Among the 118 bootstrap replicates that select Emax, suppose the \(EC_{50}\) estimates have a median of 1.35 mg/L and a 2.5th-to-97.5th percentile range of 0.62 to 3.10 mg/L. The range is conditional on Emax being selected, and its width indicates substantial instability in this parameter under the specified search. The bootstrap does not by itself prove that 1.20 mg/L is biased, nor does that conditional range have guaranteed post-selection coverage. The other models do not share the same \(EC_{50}\) interpretation, so pooling their parameter values into one distribution would be inappropriate.
Estimating Apparent-Fit Optimism
A full-search bootstrap can also estimate how much the selected procedure's apparent predictive score is inflated. For a score where larger is better, one common optimism calculation compares the selected procedure's score on the bootstrap sample used for fitting with its score when evaluated on the original data. For replicate \(b\), define \(O_b=S_{\text{train},b}-S_{\text{eval},b}\). The mean \(\bar{O}\) estimates optimism, and an optimism-corrected apparent score is \(S_{\text{apparent}}-\bar{O}\). The score, evaluation rule, and whether the selected model or the entire fitting procedure is evaluated must be specified. This calculation is an approximation, not a substitute for independent validation.
Worked Example: Correcting a Selected Predictive Score
For a hypothetical data set, a search selects one exposure-response model and gives an apparent mean predictive log score of 1.90 per subject, with higher values better. In 200 subject-level bootstrap replicates, the mean score of the selected procedure on the bootstrap training samples is 1.82 per subject, and its mean score when evaluated on the original data is 1.57 per subject. Treat these reported means as illustrative outputs from the same prespecified scoring procedure.
The estimated optimism is
Subtracting this from the apparent score gives
Thus, the selected model's apparent score is estimated to be optimistic by 0.25 units per subject under this bootstrap procedure. The corrected score is an estimate, not an observed independent-test result. Evaluation on the original data is not fully independent of the bootstrap fits, so external or properly nested validation is preferable when feasible.
Diagnostics, Alternatives, and Limits
A selection bootstrap is informative only if the computational results are trustworthy. Report the number of attempted and successful replicates, convergence rates by candidate, boundary estimates, and whether failures are concentrated in particular resamples or models. If some fits fail, do not quietly count those replicates as though they selected another model. Predefine how failures are handled, investigate their causes, and report how the conclusion changes under reasonable alternatives. A high failure rate can indicate weak identification, sparse information, or an over-complex candidate set rather than a mere software inconvenience.
The bootstrap's reliability also depends on the sampling scheme representing the target data-generating process. Subject resampling assumes subjects are reasonably exchangeable draws from the target population; stratified, clustered, or otherwise designed samples may require a resampling scheme that preserves that design. With few subjects, the empirical distribution may poorly represent population variability. Non-smooth selection rules, parameters near boundaries, and poorly identified nonlinear models can make ordinary bootstrap distributions unstable. Selection frequencies should therefore be interpreted as sensitivity to resampling the observed subjects, not as a complete account of uncertainty from model misspecification or an incomplete candidate set.
| Approach | What it addresses | Important limitation |
|---|---|---|
| Multiplicity adjustment for tests | Controls a specified family-wise error rate or false-discovery criterion across tests. | Does not by itself estimate predictive optimism or characterize the distribution of selected PK/PD parameters. |
| Full-search bootstrap | Shows how the specified model-selection procedure varies across resampled data sets; can estimate optimism. | Relies on the resampling scheme and empirical data representing the target process; does not guarantee post-selection coverage. |
| Nested cross-validation or an independent validation set | Evaluates performance while separating model selection from evaluation more directly. | Can be data-intensive and may have high variability, especially with few subjects or complex models. |
| Model averaging or predictive stacking | Combines predictions across candidates rather than relying solely on one selected model. | Does not eliminate candidate-set, extrapolation, or parameter-interpretation problems. |
These approaches answer related but distinct questions. Multiplicity control concerns error rates for a defined family of tests; model averaging and predictive stacking combine candidate predictions; validation assesses performance beyond the data used for selection. A bootstrap of the entire search focuses on how the specified selection-and-fitting procedure behaves under resampling. It should not be presented as a universal correction for every inferential consequence of selection.
A common mistake is to bootstrap only the winning model and then treat the resulting parameter spread as if the model had been fixed in advance. That procedure can describe resampling variability conditional on the chosen structure, but it omits the possibility that another candidate would win. The correction is to rerun the complete prespecified search in every replicate, report both model-selection frequencies and model-specific parameter summaries, and state clearly which summaries are conditional on selection.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Under independence, what is the probability of at least one false rejection when five null tests are each conducted at \(\alpha=0.01\)?
- Why should repeated PK/PD records generally be resampled by subject rather than by individual record?
- A candidate model is selected in 72 of 240 successful bootstrap replicates. What is its selection frequency, and what does that frequency mean?
- A bootstrap reports a parameter interval using only replicates in which one model was selected. What is the conditioning, and what should not be claimed automatically?
- In a higher-is-better score analysis, the mean bootstrap training score is 2.10 and the mean evaluation score is 1.74. What is the estimated optimism?
Solutions. (1) \(1-0.99^5=1-0.9509900499=0.0490099501\), approximately 0.04901, or 4.90%. (2) Subject-level resampling keeps each person's repeated concentration, response, dose, and covariate records together, preserving their within-subject dependence. (3) \(72/240=0.300\), or 30.0%. It is the fraction of successful resamples in which the prespecified search selected that model, not the probability that the model is true. (4) The distribution is conditional on that model winning the search. Its interval does not automatically have nominal post-selection coverage. (5) The estimated optimism is \(2.10-1.74=0.36\) score units per evaluation unit.