Worked Example: Comparing Search Strategies by Simulation

Why Compare Search Strategies by Simulation?

A covariate search procedure is judged not only by whether it finds a statistically supported model in one dataset, but also by how often it identifies the covariates that truly affect the parameter of interest. In real data, that true model is unknown. Simulation makes a controlled comparison possible: specify a population PK model with known covariate effects, repeatedly generate datasets from it, apply each search strategy, and measure what each method selects.

The result is not a universal ranking. It describes performance under the particular population model, sample size, sampling schedule, covariate distribution, residual variability, and selection rules simulated. A method that performs well when covariates are independent and effects are strong may behave differently with correlated covariates, sparse individual information, or weaker effects.

These definitions separate three questions that can otherwise be conflated. Exact recovery summarizes whether the whole selected model is correct. False inclusions measure unnecessary effects, while missed effects measure failure to retain true effects. A procedure can have a low false-inclusion count but miss several real effects, or detect most real effects while also selecting noise.

A Hypothetical Population PK Simulation

Consider a hypothetical one-compartment intravenous bolus model. For individual \(i\), dose \(D\), and sampling time \(t_{ij}\), concentration is generated as

Clearance and volume are generated as

Here \(CL_i\) is in L/h, \(V_i\) is in L, dose is in mg, time is in hours, and concentration is in mg/L. The six covariates \(x_{i1},\ldots,x_{i6}\) are standardized, independent variables with mean zero and standard deviation one. Only \(x_{i1}\) and \(x_{i2}\) affect clearance in the generating model; \(x_{i3},\ldots,x_{i6}\) are null candidates. The random effects are independent with \(\eta_i\sim N(0,0.25^2)\) and \(\zeta_i\sim N(0,0.20^2)\); the residual errors are independent with \(\epsilon_{ij}\sim N(0,0.15^2)\). Thus the residual error is lognormal on the concentration scale.

Suppose each replicate contains 240 individuals, each observed at 0.5, 1, 2, 4, 8, and 12 hours after dosing. Every candidate model uses the same structural model, random-effects structure, residual model, and fitting procedure. Candidate covariate effects are tested on clearance, not volume. This common setup is essential: otherwise differences attributed to search strategy could instead arise from different structural assumptions or data handling.

Performance Measures and Their Denominators

For each simulated dataset, record which of the six effects the procedure selects. Across replicates, count exact recoveries, false inclusions, and missed effects. The denominators matter: with 100 replicates, there are 100 opportunities for each true effect to be missed, hence 200 true-effect opportunities; there are 400 null-effect opportunities because four null covariates are considered in each replicate.

A rate based on opportunities is not the same as the proportion of replicates containing at least one error. Both can be useful, but they answer different questions. State which one is reported. Exact recovery is the number of replicates selecting precisely \(x_1\) and \(x_2\), divided by the number of replicates.

Comparing Methods on the Same Replicates

When procedures are applied to the same simulated datasets, their outcomes are paired. A direct comparison of exact recovery should preserve that pairing: each replicate provides a matched result for forward and backward selection. The difference in marginal recovery percentages alone does not show how often the two methods agree or which replicates account for their difference.

Why the Simulation Results Need Diagnostics

A selection summary is only meaningful if the simulation and fitting pipeline behave as intended. Check the generating code against the stated model, including the signs and magnitudes of the two nonzero effects and the zero effects for the remaining candidates. Check that units, sampling times, covariate standardization, and error distributions are correct. A coding error that gives a null candidate a nonzero coefficient invalidates the labels used to count false inclusions.

Also inspect fit failures by method and replicate. Record convergence status, boundary estimates, implausible parameter estimates, and whether likelihood comparisons use comparable converged fits. Quietly dropping difficult replicates can bias the comparison, especially if one selection procedure fails more often. Report how failures are handled and how many occurred. A simulation with many failed fits does not support a clean ranking of selection performance.

Monte Carlo variation is another diagnostic. Even with a fixed generating scenario, a finite number of replicates gives a noisy estimate of a method’s recovery rate. For a proportion \(\hat p\) estimated from \(R\) independent replicates, its approximate Monte Carlo standard error is \(\sqrt{\hat p(1-\hat p)/R}\). At 72% recovery with \(R=100\), this is \(\sqrt{0.72(0.28)/100}=0.0449\), or about 4.49 percentage points. That amount of simulation error cautions against treating a small difference in estimated recovery as decisive. For comparisons on shared replicates, paired uncertainty should use the replicate-level paired outcomes rather than treating the two proportions as independent.

What the Comparison Can and Cannot Establish

Forward and backward procedures use sequential testing rules. Their result can depend on entry or removal thresholds, ordering among eligible effects, and the effect of earlier selections on later tests. BIC subset selection instead searches all candidate subsets and trades fit against a penalty for model size. This exhaustive approach is feasible here because six candidates produce only 64 subsets; it may become computationally impractical with many candidates. The methods therefore differ in both search path and decision criterion.

The simulation supports a conditional statement: under the specified model and design, the recorded performance estimates were as shown. It does not prove that the generating model is scientifically correct, that BIC is generally superior, or that a selected covariate causes a change in clearance. Independence among the simulated covariates is especially consequential: correlated candidates can make effects harder to distinguish, and the performance comparison may change. Other scenarios can vary covariate correlation, effect size, sample size, residual error, sampling schedule, random-effects variability, and structural-model correctness.

A useful extension is to repeat the comparison across a planned set of scenarios, changing one or more design features at a time. The earlier material on simulation-based evaluation of covariate search procedures, power analysis, bootstrap stability, and cross-validation addresses related evaluation questions. Here, the central lesson is to make the generating truth known, define the decision rules in advance, use common replicates, and report multiple error measures with explicit denominators.

Common Pitfalls

• Using only exact recovery. Two procedures can have similar exact-recovery proportions but different balances of false inclusions and missed effects. Report each performance measure relevant to the decision.
• Confusing opportunities with datasets. A false-inclusion rate based on 400 null-covariate opportunities is not the same as the percentage of 100 datasets with at least one false inclusion. Name the denominator.
• Ignoring pairing. The procedures were applied to the same replicates. Their concordance and discordance should be counted on those matched datasets, not inferred from marginal percentages.
• Generalizing beyond the scenario. A result from independent standardized covariates and a particular sampling design does not settle performance under collinearity or sparse sampling.
• Counting failed fits as ordinary selections. Convergence failures require explicit reporting and a prespecified handling rule; they are not automatically correct model selections.