Predictive Stacking for PK/PD Model Ensembles

Why Stack PK/PD Predictions?

Competing PK/PD models can make different predictions for the same subject, dose, and time point. One model may describe an exposure-response curve well, while another captures a delayed effect or a different residual structure. Selecting one model discards the others; averaging can retain predictive contributions from several. Predictive stacking chooses the weights in that average to optimize out-of-sample predictive density, rather than to represent the probability that each candidate model is true.

The prediction target matters. In repeated-measures PK/PD data, a useful target is often the response trajectory of a new subject, conditional on that subject’s known design information, such as dose and sampling times. For that target, validation must hold out whole subjects: holding out one record while using the same subject’s other records to fit the model would answer a different question.

Each weighted sum inside the logarithm is the ensemble’s predictive density for one held-out subject. Thus the objective is the sum of subject-level log predictive densities. The density must refer to the same response data, units, and prediction target across candidate models. If a response transformation is used, densities must be put on a common scale, including the appropriate transformation adjustment.

Why the Stacking Objective Is Concave

For two models, write \(w\) for the weight on model 1 and \(1-w\) for the weight on model 2, where \(0\leq w\leq1\). The objective is $$ L(w)=\sum_{i=1}^{n}\log\left(wq_{1i}+(1-w)q_{2i}\right). $$ Assume the predictive densities are positive and finite. Differentiating gives $$ L'(w)=\sum_{i=1}^{n}\frac{q_{1i}-q_{2i}}{wq_{1i}+(1-w)q_{2i}}, $$ and a second differentiation gives $$ L''(w)=-\sum_{i=1}^{n}\frac{(q_{1i}-q_{2i})^2}{\left(wq_{1i}+(1-w)q_{2i}\right)^2}\leq 0. $$

This establishes concavity: a maximum can be found reliably on the feasible interval. If at least one subject has different predictive densities under the two models, the second derivative is strictly negative throughout the interval. The objective is then strictly concave and has a unique maximum on \([0,1]\). That maximum may be at a boundary, in which case there need not be an interior stationary point. With more than two models, the objective is concave in the weight vector on the simplex; numerical constrained optimization is commonly used.

Stacking Weights and Posterior Model Probabilities

Posterior model probabilities and stacking weights both form nonnegative weights summing to one, but they answer different questions. A posterior model probability is \(P(M_m\mid D)\), calculated from a prior probability on each model and its marginal likelihood. Under Bayesian model averaging, those probabilities weight model-specific predictive distributions. Stacking instead selects weights to maximize the sum of leave-one-subject-out log predictive densities. Its weights are not posterior probabilities and do not state the probability that a candidate model is true.

This distinction is useful when no candidate model is expected to represent every feature of the data. A model with a relatively low posterior probability may still contribute useful predictions for some held-out subjects. Conversely, a model with high posterior probability under a particular set of priors need not receive the largest stacking weight. The posterior model probabilities are sensitive to prior model odds and marginal likelihoods; stacking directly targets the chosen predictive score. Neither approach is automatically superior: the choice depends on whether the scientific target is model uncertainty under a specified Bayesian model set or predictive performance for a specified validation target.

Calculating an Ensemble Prediction and Its Uncertainty

For a future response \(Y\), let model \(m\) have predictive mean \(\mu_m\) and predictive variance \(\sigma_m^2\), conditional on the same covariates and prediction target. The stacked predictive distribution is the mixture \(\sum_m w_m p_m(Y)\). Its mean is \(\sum_m w_m\mu_m\), and its variance includes both within-model predictive variance and between-model variation in predictive means: $$ \operatorname{Var}(Y)=\sum_{m=1}^{M}w_m\left[\sigma_m^2+(\mu_m-\mu_{\mathrm{ens}})^2\right], \qquad \mu_{\mathrm{ens}}=\sum_{m=1}^{M}w_m\mu_m. $$ This follows by taking the weighted average of each component’s second moment, \(\sigma_m^2+\mu_m^2\), and subtracting \(\mu_{\mathrm{ens}}^2\). The resulting mixture need not itself be normal, even when every component predictive distribution is normal.

Practical Conditions, Diagnostics, and Limitations

A stacking analysis is interpretable only when its validation design matches the intended prediction. For a new-subject target, every record from a held-out subject must be excluded from fitting, including information used to estimate that subject’s random effect. The predictive density should then integrate over the random-effect distribution for a new subject. If the actual target is a future observation from a subject whose earlier observations are available, a conditional prediction may be appropriate instead; that is a different target and should use a matching validation design.

Before trusting the fitted weights, check that each fold converged, the predictive densities are finite and computed on a common scale, and the weights are nonnegative and sum to one. Inspect the subject-level log-score contributions and the fitted weights. A boundary solution can be a valid optimum, but it can also indicate that one candidate dominates under this scoring target or that the candidate set lacks useful diversity. Large changes in weights after small changes to the candidate set or held-out subjects signal instability. These checks diagnose fragility; they do not by themselves establish good calibration or reliable extrapolation.

The scoring unit also matters. Giving each subject one joint trajectory density gives each subject one log-score contribution, but a long, densely sampled trajectory can still have a very different joint density scale from a short trajectory. Decide whether the target is the full trajectory, particular time points, or another defined prediction unit before constructing the score. The log score is sensitive to assigning very small density to an observed outcome, so a single poorly predicted subject can have substantial influence.

Stacking depends on the candidate set and on the observed validation cases. It cannot rescue an ensemble if all candidate models fail in the same way, and optimized weights do not guarantee performance outside the covariate, dose, or time ranges represented by the validation data. A posterior model probability offers a coherent Bayesian summary conditional on its models and priors; stacking offers a data-driven predictive combination for a specified scoring target. A single selected model is simpler to interpret and deploy, but may discard predictive information. These are alternatives with different aims, not interchangeable estimates of one quantity.