Tutorials › Pharmacokinetic-Pharmacodynamic (PK/PD) Modeling › Subject-Level Cross-Validation for Repeated PK/PD Data

Model Selection and Averaging · Tutorial 479 of 598

Subject-Level Cross-Validation for Repeated PK/PD Data

Use leave-one-subject-out validation to assess how well competing PK/PD models predict the outcomes of patients not used for fitting.

Advanced 12 min read

What You'll Learn

  • Define the new-subject prediction target for repeated PK/PD data
  • Construct leave-one-subject-out training and validation folds
  • Distinguish population predictions from predictions informed by a held-out subject’s outcomes
  • Calculate and interpret joint predictive log scores
  • Compare subject-level validation with observation-level splitting and information criteria
  • Identify leakage, dependence, and unequal follow-up as sources of misleading scores

Why the Validation Unit Matters

Repeated PK/PD measurements from one patient are not independent pieces of evidence about a new patient. They share that patient’s dose and covariate history, and may share latent characteristics represented by subject-specific random effects. If the goal is to predict outcomes for a patient who was not used to fit the model, validation must keep all measurements from each patient together.

Subject-level cross-validation evaluates predictions for held-out subjects. In leave-one-subject-out (LOSO) validation, each subject is held out once, the model is refitted using all other subjects, and predictions are generated for the held-out subject. This targets performance on new subjects drawn from a population like the study population, under the specified prediction conditions. It is different from asking how well a model predicts another measurement from a patient whose earlier measurements are already available.

Definition: Let \(y_i=(y_{i1},\ldots,y_{in_i})\) denote all observed outcomes for subject \(i\), with \(n_i\) measurements, and let \(D_{-i}\) be the data from every subject except \(i\). Subject-level LOSO validation fits each candidate model using \(D_{-i}\), then evaluates predictions for the complete vector \(y_i\) without using any of subject \(i\)’s outcomes to fit the model or estimate that subject’s random effects.

The held-out subject’s covariates and dosing information may be used if they would be available when making the intended prediction. Outcomes from that subject must not be used to tune the model, select covariates, estimate subject-specific effects, or choose among candidate structures within that fold. A different prediction task may legitimately condition on a patient’s earlier outcomes, but that requires a validation design that represents the information available at the forecast time.

Constructing LOSO Predictions

Suppose a candidate nonlinear mixed-effects model has population parameters \(\theta\), subject random effects \(b_i\), and an observation model for \(y_i\). The predictive distribution for a new subject’s complete outcome vector averages over the random effects that are unknown for that subject. In a Bayesian analysis it also averages over parameter uncertainty given the training data:

$$ p(y_i\mid D_{-i},M) = \int p(y_i\mid b_i,\theta,M)\, p(b_i\mid\theta,M)\, p(\theta\mid D_{-i},M)\,db_i\,d\theta. $$

Here \(M\) denotes the candidate model. In a frequentist analysis, a common practical approximation estimates \(\theta\) from \(D_{-i}\), then integrates over the fitted population distribution of \(b_i\):

$$ p(y_i\mid D_{-i},M) \approx \int p(y_i\mid b_i,\hat{\theta}_{-i},M)\, p(b_i\mid\hat{\theta}_{-i},M)\,db_i. $$

This is a population-level prediction: it does not insert an empirical Bayes estimate of \(b_i\) calculated from the held-out subject’s outcomes. For a nonlinear model, the integral may require numerical quadrature or simulation. A plug-in estimate of population parameters does not account fully for their estimation uncertainty; where that uncertainty matters, it can be included using a suitable parameter-sampling procedure.

1
Define the prediction target.
Specify whether the aim is to predict a new subject’s full trajectory, a particular endpoint, or future outcomes conditional on earlier observations. State which covariates and dosing information are available at prediction time.
2
Hold out one complete subject.
Remove that subject’s entire outcome history. Keep all repeated measurements for the remaining subjects in the training set.
3
Refit and predict.
Fit the candidate model on the training subjects, then obtain a predictive distribution or prediction for the held-out subject using only information allowed by the target.
4
Score and aggregate.
Calculate the prespecified score for the held-out subject. Repeat for every subject and summarize the subject-level scores using a consistent rule.

Every data-dependent step belongs inside the training fold. This includes estimating preprocessing quantities, selecting covariates from outcomes, or choosing a model based on fit. If candidate structures are fixed in advance, each one is still fitted anew in each fold. Convergence failures should be recorded and investigated rather than silently dropped, since excluding difficult folds can make a model look artificially favorable.

Why Observation-Level Splitting Leaks Subject Information

If some observations from a subject enter the training set while other observations from that same subject are treated as test data, the model has information about the held-out subject. This is especially consequential when the model includes random effects: training outcomes help estimate that subject’s latent effect, which then improves predictions for the test observations. Such a split can answer a useful within-subject forecasting question, but it does not estimate prediction for a new patient with no outcome history.

A simple normal random-intercept model shows the mechanism. Suppose \(b_i\sim N(0,4)\), and one measurement is \(y_{i1}=b_i+\epsilon_{i1}\), where \(\epsilon_{i1}\sim N(0,1)\) independently. After observing a training value \(y_{i1}=10\), the conditional mean of the subject effect is:

$$ E(b_i\mid y_{i1}=10) = \frac{4}{4+1}(10) = 8. $$

A subsequent measurement from that same subject is therefore predicted using a subject-effect estimate near 8, rather than the new-subject population mean of 0. The information in the first measurement is legitimate for forecasting within that patient, but it makes an observation-split score optimistic for the new-patient target.

Worked Example: Assigning Folds Without Splitting Patients

A hypothetical PK/PD dataset contains 12 subjects, each with four concentration or response measurements. For a three-fold comparison, assign subjects 1–4 to fold A, 5–8 to fold B, and 9–12 to fold C.

When fold A is tested, all 16 observations from subjects 1–4 are withheld together. The model is fitted using all 32 observations from subjects 5–12, and predictions are made for each of the four held-out subjects. Folds B and C are handled in the same way. Every subject is tested once, and no subject contributes outcomes to both the training and test data in a fold.

A random split of 48 individual measurements could, for example, put three measurements from subject 1 in training and the fourth in testing. That test measurement would then be predicted with information about subject 1’s response history. It would not represent a prediction for a previously unseen patient.

Choosing and Interpreting Predictive Scores

A score must match the intended prediction target. For a complete repeated-measures trajectory, a natural choice is the joint log predictive density, \(\log p(y_i\mid D_{-i},M)\). Larger values indicate that the model assigned greater predictive density to the observed held-out trajectory. Equivalently, one can report the negative log predictive density, for which smaller is better. The joint score should reflect the model’s dependence structure; treating correlated measurements as independent changes the score and can misrepresent predictive performance.

For a target focused on point predictions, one may instead compare subject-level errors, such as the mean squared error over a prespecified set of measurements. If each subject is to count equally, first summarize error within each subject and then average across subjects. Pooling every observation gives more weight to subjects with more measurements. Neither weighting is universally correct: the choice should represent the population and decision the analysis is intended to support.

Worked Example: Calculating a Joint Predictive Log Score

Consider a hypothetical two-measurement trajectory for one held-out subject, \(y_i=(2,3)\), under a random-intercept model with \(b_i\sim N(0,4)\) and independent measurement errors with variance 1. After integrating over the new subject’s random intercept, the predictive covariance matrix is:

$$ \Sigma = \begin{pmatrix} 5 & 4\\ 4 & 5 \end{pmatrix}, \qquad |\Sigma|=9, \qquad \Sigma^{-1} = \frac{1}{9} \begin{pmatrix} 5 & -4\\ -4 & 5 \end{pmatrix}. $$

The predictive mean is \(\mu=(0,0)\). For a bivariate normal density, the log predictive density is \(-\log(2\pi)-\frac{1}{2}\log|\Sigma|-\frac{1}{2}(y_i-\mu)^\mathsf{T}\Sigma^{-1}(y_i-\mu)\). Substituting the values gives:

$$ \begin{aligned} (y_i-\mu)^\mathsf{T}\Sigma^{-1}(y_i-\mu) &= \frac{5(2^2)-8(2)(3)+5(3^2)}{9} =\frac{17}{9} \approx 1.889,\\ \log p(y_i\mid M) &= -\log(2\pi)-\frac{1}{2}\log(9)-\frac{1}{2}\left(\frac{17}{9}\right)\\ &\approx -1.838-1.099-0.944 =-3.881. \end{aligned} $$

Thus this model assigns log predictive density \(-3.881\) to the held-out two-measurement vector. The value is not a probability that the model is correct. It is a density-based score on the scale of the observed outcomes; comparisons are meaningful when models score the same subjects and the same outcome vectors under the same scoring convention.

If the two observations were incorrectly treated as independent with marginal variance 5, the log density would be \(-\log(2\pi)-\log(5)-(2^2+3^2)/(2\cdot5)\approx-4.747\). This difference illustrates that a joint score depends on the assumed covariance, not merely on separate marginal predictions.

Worked Example: Comparing Two Models Across Subjects

Suppose two hypothetical models are evaluated using the same five held-out subjects. Define each difference as the Model A log score minus the Model B log score. The subject-level differences are \(0.4,-0.2,0.3,0.1,-0.1\) nats per subject.

The mean difference is:

$$ \overline{\Delta} = \frac{0.4-0.2+0.3+0.1-0.1}{5} = 0.10\text{ nats per subject}. $$

The sample standard deviation of the five differences is \(s_\Delta=\sqrt{0.26/4}\approx0.255\) nats per subject, so the descriptive standard error is \(s_\Delta/\sqrt{5}\approx0.114\) nats per subject. Model A has the higher average score in this example, but the difference is small relative to the observed variation across subjects. These summaries do not establish that Model A is generally superior; with five subjects, the comparison is imprecise, and the folds’ fitted models also share much of their training data.

Assumptions, Diagnostics, and Alternatives

LOSO results are most informative when subjects are independent sampling units, the held-out subjects represent the intended prediction population, and the candidate models are fitted successfully in each fold. Within-subject dependence need not be absent; it should be represented appropriately by the model and retained in the held-out subject’s joint score. The prediction target and scoring rule should be fixed before comparing candidates. If the data contain site, family, or other higher-level clusters that create dependence across subjects, holding out single subjects may still leak cluster-level information; the held-out unit may need to be the relevant cluster.

Inspect fold-level scores and predictions, not only their average. A few subjects may dominate a joint log score, particularly if they have long follow-up or extreme trajectories. Check whether poor scores coincide with unusual covariates, sparse sampling, unexpected dosing histories, convergence problems, or implausible predictive distributions. For each fold, verify that the predicted distribution is finite and that the observed outcomes are on the scale and within the support specified by the model. Report how many subjects and measurements contribute, and explain how subjects with different numbers of observations are weighted.

ApproachWhat it evaluatesKey distinction
Subject-level LOSOPrediction for a subject omitted from fittingHolds out the complete subject trajectory
Observation-level splittingPrediction of measurements when other measurements from the same subject are availableCan use subject-specific information; not a new-subject assessment
AICc or BICPenalized in-sample likelihood for candidate modelsDoes not directly score predictions for held-out subjects

Information criteria can be useful for comparing candidate structures, but they answer a different question from cross-validated prediction. LOSO is also computationally demanding because the model must be refitted for each held-out subject, and scores can be unstable with few subjects or influential trajectories. It estimates predictive performance only for the population, observation process, and prediction horizon represented by the validation design; it does not guarantee performance under a new dose range, site, or sampling scheme.

Common Pitfall

Key takeaway: Do not split repeated measurements at random and then describe the resulting score as performance for new patients. A test observation from a subject represented in training can benefit from that subject’s estimated random effect. Hold out all outcomes from a subject, integrate over the random-effect distribution for new-subject prediction, and interpret scores relative to a clearly stated target.

Another common error is to compare totals formed from different numbers of observations. A subject with a long trajectory can contribute a larger-magnitude joint log score than a subject with only a few measurements. This is not automatically a flaw, but it makes the aggregation rule part of the estimand. Decide whether the target is predictive quality per patient, per measurement, or for the complete trajectory, and report the choice.

Check Your Understanding

Answer each question, then compare with the brief solutions.

  1. For a new-patient prediction target, should five observations from one subject be split between training and testing? Why?
  2. In a random-effects model, what must be done with the random effect for a held-out subject whose outcomes are unavailable for fitting?
  3. Two models have mean LOSO log scores of \(-3.2\) and \(-3.7\) on the same subjects. Which has the higher average score, and what does that comparison mean?
  4. Why can a pooled observation-level error give a different comparison from an average of subject-level errors?

Solutions. (1) No. Hold out all five outcomes together; otherwise the training data contain information about that patient and the test does not represent a new patient. (2) Integrate over the fitted population distribution of random effects, rather than estimating that subject’s random effect from held-out outcomes. (3) The model scoring \(-3.2\) has the higher average log predictive density and assigned greater predictive density to the observed held-out outcomes on average; this alone does not prove general superiority. (4) Pooled error weights subjects with more observations more heavily, while averaging subject-level errors can give each patient equal weight. These summarize different targets when follow-up counts vary.

Clinical Biostats