Cross-Validation for Covariate Model Selection in Population PK

Why Validate Covariate Models by Subject?

A covariate model can improve the fit to the data used to build it without improving predictions for new patients. Cross-validation addresses this distinction by repeatedly fitting candidate models on one set of subjects and assessing their predictions on different subjects. In population pharmacokinetics (PK), the split must respect the sampling structure: records from the same patient share that patient's latent random effects and cannot be treated as independent validation units.

The target considered here is prediction for a new subject drawn from the same target population, using that subject's covariates and sampling design but not using their concentration observations to fit the model. This is a different prediction task from forecasting later concentrations for a patient whose earlier concentrations are already available. The latter may use those earlier records, but only if the validation design explicitly defines that information as available at prediction time.

What Is Fitted and What Is Scored?

A typical population PK model describes an observation \(Y_{ij}\) for subject \(i\) at occasion or time \(j\) using a structural model, subject-level random effects, and a residual-error model. Candidate covariate models may differ in which covariates they include or in the functions used to represent covariate effects. In each training fold, all parameters that belong to the candidate fitting procedure are estimated using training subjects only.

For a held-out subject with covariates \(x_i\) and concentration records \(y_i\), a suitable new-subject predictive density is the population density obtained by integrating over that subject's unobserved random effects. If \(\theta\) denotes population parameters, \(\eta_i\) the subject random effects, and \(\Omega\) their distributional parameters, a plug-in predictive density has the form:

Here \(M\) denotes the candidate model and the estimates are obtained from the training subjects in that fold. The integral matters: a new subject's random effect is unknown before observing that subject's concentrations. Substituting an empirical Bayes estimate obtained from the validation concentrations instead gives the scoring procedure information it would not have for a genuinely new patient.

One common score is the subject-level log predictive density, \(\ell_i=\log p(y_i\mid x_i,\widehat{\theta}_{\mathrm{train}},\widehat{\Omega}_{\mathrm{train}},M)\). Larger values indicate that the model assigned greater predictive density to the observed validation records. Another option is a prediction-error measure, such as squared error on a prespecified concentration scale. The prediction and the unit of averaging must match the intended use. Log predictive densities based on different numbers of records are totals over different amounts of information; averaging total scores equally over subjects and dividing the total score by the number of records answer different questions.

A Subject-Level Cross-Validation Procedure

In \(K\)-fold cross-validation, subjects are partitioned into \(K\) groups. Each group is held out once, while the other \(K-1\) groups form the training set. A complete comparison uses the same folds, scoring rule, and prediction target for every candidate model.

Why the Held-Out Score Estimates Prediction

The justification is an averaging argument. Suppose subjects are independent draws from the target population and the folds are formed without using their outcomes. For a fixed fitted training model, each held-out subject represents a draw from the same target distribution as a future subject. Its score therefore provides information about the model's predictive performance on that distribution. Averaging held-out scores across subjects estimates the corresponding average performance.

There is an important qualification: the model fitted in each fold uses fewer subjects than a final model fitted to the full data. With \(K\) approximately equal folds, each fold-specific fit uses roughly \((K-1)n/K\) subjects. The cross-validation result therefore reflects the specified fitting procedure at approximately that training-set size; it is not an exact guarantee of how a full-data fit will perform.

The justification also depends on the split respecting dependence. If the same patient's records appear in both training and validation sets, training data reveal information about that patient's random effects. The validation records then no longer represent prediction for an unseen subject. Similarly, preprocessing or covariate screening performed once on all outcomes allows held-out information to influence the fit, even if the final parameter estimation is repeated within folds.

Worked Example: Detecting Leakage From Validation Concentrations

Consider a simplified hypothetical random-intercept model for one concentration per subject: \(Y_i=\mu+\eta_i+\epsilon_i\), with \(\mu=0\), \(\eta_i\sim N(0,1)\), and \(\epsilon_i\sim N(0,1)\), independently. For a new subject, before observing their concentration, the marginal distribution is \(Y_i\sim N(0,2)\). Suppose the held-out value is \(y_i=2\).

A valid new-subject log predictive density under this model is:

Now consider the incorrect procedure of using the held-out value to estimate the subject's random effect before scoring that same value. The conditional mean of \(\eta_i\) given \(Y_i=2\) is \(1\), because the two independent normal components have equal variance. If this estimate is plugged in and the observation is scored using the residual model \(N(\eta_i,1)\), the score becomes:

This score is higher, but the apparent improvement is artificial: the concentration being scored was used to choose the random-effect estimate. A correct new-subject assessment integrates over the random-effect distribution rather than treating a validation-derived estimate as known. This example is deliberately simple; in a nonlinear PK model, the same leakage can arise when validation concentrations are used to calculate empirical Bayes estimates.

Diagnostics, Alternatives, and Limitations

A cross-validation result is only as credible as its split and fitting pipeline. Before interpreting model rankings, check that every subject appears in validation exactly once for ordinary \(K\)-fold evaluation, that no subject crosses the training-validation boundary, and that candidate models use identical folds. Inspect fit convergence, parameter estimates, and predictive-score distributions across folds and subjects. A result driven by one unusually influential subject, a failed fold fit, or a small number of extreme concentrations deserves investigation rather than automatic acceptance.

• Information leakage: full-data screening, outcome-based transformation choices, imputation, or scaling can transmit validation information into training. Re-estimate data-dependent steps within each training fold.
• Unrepresentative folds: small samples or rare covariate groups can leave a fold with little support for a covariate effect. Examine covariate distributions and avoid splitting subjects in ways that make a fold's training data unable to estimate the candidate model.
• Model instability: nonlinear mixed-effects fits may fail or converge to different solutions across folds. Report how failures are handled, and do not silently omit difficult folds.
• Score mismatch: a score on observed concentrations reflects the chosen residual model and observation schedule. It does not by itself measure dosing decisions, clinical benefit, or performance in a population with a different sampling design.

Cross-validation complements, rather than replaces, other approaches. Likelihood-based comparisons evaluate fit to the analyzed data under their assumptions; information criteria adjust fit using a complexity penalty, while cross-validation directly evaluates held-out prediction. Bootstrap stability summarizes how often selection changes under resampling, whereas subject-level cross-validation compares predictive scores on held-out subjects. A single train-validation split is simpler but can be highly dependent on the chosen partition. External validation is often a stronger assessment of transport to a distinct population when suitable data are available.

If many covariate models are compared and the one with the best cross-validation score is selected, that best score is itself subject to selection optimism: random favorable variation helped determine which model won. The score is useful for comparing candidates under a common procedure, but should not be presented as an independent, unbiased estimate of the selected model's final performance. A separate untouched test set or nested cross-validation can provide a more independent assessment; both require enough subjects to support the additional split or fitting steps. Cross-validation also cannot repair misspecification shared by every candidate, correct biased or unrepresentative data, or establish that an observed covariate association is causal.

Common Mistake

Wrong reasoning: “I split concentration records at random, fitted the model on most records, and used empirical Bayes estimates from the held-out records to calculate validation predictions. The validation score is independent because those records were not used to estimate the population parameters.”

Correction: records from one subject are dependent through that subject's random effects, and validation-derived empirical Bayes estimates use the outcomes being scored. For new-subject prediction, assign whole subjects to folds and integrate over the random-effect distribution when scoring their concentrations. Keep all data-driven model-building steps inside each training fold as well.