Bayesian Model Averaging for NLME Predictions

Why Average Predictive Distributions?

External prediction often depends on choices about the structural model, covariate effects, or random-effects specification. Selecting one model and treating it as certainly correct can hide uncertainty that matters for an external subject’s prediction. Bayesian model averaging (BMA) represents uncertainty across a specified set of candidate models by combining their predictive distributions.

The object being averaged is the distribution of the future external outcome, not the model parameters. This builds on the earlier material on averaging predictive distributions: each candidate NLME model first supplies a predictive distribution for the same target, and the model probabilities determine how those distributions are combined. The resulting distribution can reflect both uncertainty within each model and differences between models.

Definition: Let $M_1,\ldots,M_K$ be candidate models and let $D$ denote the data used to fit and compare them. The posterior model probability for $M_k$ is $w_k=P(M_k\mid D)$, where $w_k\geq 0$ and $\sum_{k=1}^K w_k=1$. For an external outcome $Y^*$ and its prediction conditions $x^*$, the Bayesian model-averaged predictive distribution is $$ p(y^*\mid x^*,D)=\sum_{k=1}^K w_k\,p(y^*\mid x^*,D,M_k). $$

Here, $x^*$ denotes the external subject’s prediction conditions, such as the relevant covariates, dose history, and prediction time. A candidate model’s predictive distribution should already account for the uncertainties required by the target: for example, posterior uncertainty in its parameters and, for a new subject, the model’s distribution for that subject’s random effects and residual error. A population-mean target and a future individual observation are different targets, so their component predictive distributions need not be the same.

For the ordinary posterior model probabilities above, the candidate models are assigned prior probabilities $P(M_k)$. Bayes’ rule gives

w_k=P(M_k\mid D) =\frac{p(D\mid M_k)P(M_k)} {\sum_{j=1}^K p(D\mid M_j)P(M_j)}.

The marginal likelihood $p(D\mid M_k)$ averages the likelihood over the parameter prior within model $M_k$. Thus, posterior model probabilities depend on both the data and the prior specification, including the model priors and parameter priors used to compute marginal likelihoods. They are not simply measures of how well each model fits a training data set at its best-fitting parameter values.

Constructing the Mixture

Before combining predictions, check that every model predicts the same outcome on the same scale, for the same external target, and conditional on the same information. The weights are based on the specified training data $D$. If an external outcome is meant to be held out for validation, it must not be used to calculate the weights or select the candidate models; doing so would compromise the external evaluation.

Once the component predictive distributions and weights are available, the mixture gives the probability of an event by weighting the corresponding probabilities under each model. For a threshold $c$, for example, the probability that the external outcome is no greater than $c$ is the weighted sum of the model-specific probabilities.

P(Y^*\leq c\mid x^*,D) =\sum_{k=1}^K w_kP(Y^*\leq c\mid x^*,D,M_k).

Worked Example: Updating Model Weights and Combining Predictions

Suppose two candidate NLME models have prior probabilities 0.50 each. Their marginal likelihoods for the training data are 0.12 and 0.03. For a hypothetical external prediction target, suppose their posterior predictive distributions have means 10.0 and 14.0 mg/L and variances 1.00 and 4.00 $(\text{mg/L})^2$, respectively. Calculate the posterior model probabilities and the mean and variance of the model-averaged prediction.

Posterior model probabilities. The unnormalised weights are $0.12(0.50)=0.060$ for $M_1$ and $0.03(0.50)=0.015$ for $M_2$. Their sum is 0.075, so

w_1=\frac{0.060}{0.075}=0.800,\qquad w_2=\frac{0.015}{0.075}=0.200.

Predictive mean. The mixture mean is the weighted average of the component means:

\mu_{\mathrm{mix}} =0.800(10.0\ \text{mg/L})+0.200(14.0\ \text{mg/L}) =10.8\ \text{mg/L}.

Predictive variance. The within-model contribution is $0.800(1.00)+0.200(4.00)=1.60\ (\text{mg/L})^2$. The model means differ from the mixture mean by $-0.8$ and $3.2$ mg/L, giving a between-model contribution of $0.800(0.8^2)+0.200(3.2^2)=2.56\ (\text{mg/L})^2$. Therefore the total predictive variance is $1.60+2.56=4.16\ (\text{mg/L})^2$, with standard deviation $\sqrt{4.16}=2.040$ mg/L. The relatively large between-model contribution arises because the candidate models predict different means.

Why Between-Model Uncertainty Appears

The mixture variance has a useful general decomposition. Write $\mu_k=E(Y^*\mid x^*,D,M_k)$ and $\sigma_k^2=\operatorname{Var}(Y^*\mid x^*,D,M_k)$. The mixture mean is $\mu_{\mathrm{mix}}=\sum_k w_k\mu_k$. Applying the law of total variance, with the model index as the conditioning variable, gives

\begin{aligned} \operatorname{Var}(Y^*\mid x^*,D) &=E\{\operatorname{Var}(Y^*\mid x^*,D,M)\mid x^*,D\}\\ &\quad+\operatorname{Var}\{E(Y^*\mid x^*,D,M)\mid x^*,D\}\\ &=\sum_{k=1}^K w_k\sigma_k^2 +\sum_{k=1}^K w_k(\mu_k-\mu_{\mathrm{mix}})^2. \end{aligned}

The first term averages the predictive variances within models. These may include parameter, random-effect, and residual uncertainty according to the prediction target. The second term measures the spread of the model-specific means around the mixture mean. It is nonnegative and is zero only when all positively weighted models have the same predictive mean. The decomposition is exact for the stated component distributions; it does not require the mixture itself to be normal.

Worked Example: A Mixture Interval Can Span Separate Modes

Consider a hypothetical predictive mixture for an external concentration: $0.800N(10,1)+0.200N(30,1)$, where each normal distribution’s second argument is its variance and concentration is measured in mg/L. Find an equal-tailed 95% predictive interval. The models’ means are far apart, so an interval based on a single normal approximation would obscure the mixture’s shape.

The mixture cumulative distribution function is $F(y)=0.800\Phi(y-10)+0.200\Phi(y-30)$, where $\Phi$ is the standard normal cumulative distribution function. At the lower 2.5% quantile, the second component’s contribution is negligible because the proposed endpoint is many standard deviations below 30 mg/L. Thus $0.800\Phi(y-10)\approx0.025$, so $\Phi(y-10)\approx0.03125$. The corresponding standard normal quantile is approximately $-1.863$, and the lower endpoint is $10-1.863(1)=8.137$ mg/L.

At the upper 97.5% quantile, the first component’s upper tail is negligible. The equation is approximately $0.800+0.200\Phi(y-30)=0.975$, so $\Phi(y-30)=0.875$. Since $\Phi^{-1}(0.875)\approx1.150$, the upper endpoint is $30+1.150(1)=31.150$ mg/L. The equal-tailed interval is therefore approximately [8.137, 31.150] mg/L. Its wide span reflects substantial uncertainty about which model is more appropriate; most of the mixture probability lies near one of the two component means rather than uniformly throughout the interval.

This is a posterior predictive interval for the specified mixture, conditional on the data, candidate models, and prior choices. Under that predictive distribution, the interval contains 95% of the probability for the future outcome. It is not a confidence interval for a model parameter.

Prediction Differences Across Several Models

The model-averaged predictive mean and variance can be calculated from component summaries when those summaries are available. However, interval endpoints generally cannot be calculated by averaging component endpoints: the mixture’s quantiles are found from its full cumulative distribution. Numerical evaluation of the weighted component CDF, or simulation from the already constructed component predictive distributions, can be used to obtain quantiles when no closed form is available.

Worked Example: Combining Three External Predictions

For one hypothetical external subject and one prediction time, three fitted NLME models have posterior model probabilities 0.50, 0.30, and 0.20. Their component predictive means are 5.0, 8.0, and 7.0 mg/L, and their predictive variances are 1.00, 1.00, and 4.00 $(\text{mg/L})^2$. Find the model-averaged mean and variance.

The weights sum to $0.50+0.30+0.20=1.00$. The predictive mean is

\mu_{\mathrm{mix}} =0.50(5.0)+0.30(8.0)+0.20(7.0) =6.30\ \text{mg/L}.

The within-model variance contribution is $0.50(1.00)+0.30(1.00)+0.20(4.00)=1.60\ (\text{mg/L})^2$. The between-model contribution, using differences from 6.30 mg/L, is

0.50(5.0-6.30)^2+0.30(8.0-6.30)^2+0.20(7.0-6.30)^2 =0.845+0.867+0.098=1.810\ (\text{mg/L})^2.

Thus the total variance is $1.60+1.81=3.41\ (\text{mg/L})^2$, and the standard deviation is $\sqrt{3.41}=1.847$ mg/L. This calculation combines the models’ predictive distributions for the same external subject and target; it does not imply that the subject’s concentration is normally distributed with this mean and variance.

Assumptions, Diagnostics, and Alternatives

BMA is conditional on the candidate model set and the prior specifications. Its probabilities describe uncertainty among those models, not the probability that one of them must be the true data-generating mechanism. If all candidates miss important structure, the model-averaged prediction can still be poorly calibrated. Likewise, posterior model probabilities can be sensitive to prior choices, especially for marginal likelihoods, so the priors and the model set should be scientifically defensible.

Useful checks include confirming that posterior model probabilities are nonnegative and sum to one; verifying that every model’s predictive distribution addresses the same outcome and external conditions; and inspecting how much each model contributes to the mixture. Large differences between model-specific predictions create substantial between-model variance and may yield a multimodal or unusually broad predictive distribution. Assess external predictive calibration and interval coverage as appropriate to the prespecified target; good agreement in a training-data model comparison alone does not establish external predictive adequacy.

The main alternatives answer somewhat different questions. Selecting a single model discards model uncertainty. Averaging parameter estimates can be incoherent when models have different parameterizations or structures, and in a nonlinear model the prediction evaluated at an averaged parameter need not equal an average of predictions. Predictive-score stacking, covered in the earlier material on stacking weights, chooses weights to optimize predictive performance rather than using posterior model probabilities. Stacking can be useful when the candidate set is viewed as an approximation rather than as a set of competing probability models, but its weights have a different interpretation. Neither approach repairs a shared failure of the candidate models or justifies using external outcomes that were supposed to remain held out.

Common pitfall: Averaging model-specific 95% interval endpoints and calling the result a 95% model-averaged interval is generally incorrect. Quantiles are nonlinear: first combine the predictive CDFs using the model weights, then find the quantiles of that mixture. Also, do not average model parameters and assume the resulting parameter value produces the model-averaged prediction.

Define the target.
Specify the external outcome and conditioning information, including whether the prediction concerns a population quantity or a future individual observation.

Establish weights.
Use the training data and stated prior probabilities to calculate posterior model probabilities; do not use held-out external outcomes.

Combine predictions.
Form the weighted predictive distribution, then calculate its summaries and quantiles from that distribution.

Assess adequacy.
Check model contributions, predictive dispersion, external calibration, and sensitivity to the candidate set and prior choices.

Check Your Understanding

Answer each question, then compare with the concise solutions below.

Two models have prior probabilities 0.60 and 0.40 and marginal likelihoods 0.10 and 0.05. What are their posterior model probabilities?
Why can a mixture predictive variance exceed the weighted average of its component variances?
Is the weighted average of two model-specific 95% interval endpoints generally a 95% interval for the mixture? Explain.
One candidate model predicts a population mean concentration and another predicts a future individual observation. Can their predictive distributions be combined directly? Why or why not?
What is the key difference in interpretation between posterior model averaging weights and predictive-score stacking weights?

Solutions. (1) The unnormalised weights are 0.060 and 0.020, which sum to 0.080, so the posterior probabilities are $0.060/0.080=0.750$ and $0.020/0.080=0.250$. (2) The total variance includes the nonnegative between-model variance of the component means, in addition to the weighted within-model variances. (3) No. Mixture quantiles are obtained from the weighted mixture CDF, not by averaging component quantiles. (4) Not directly: the models must predict the same outcome on the same scale for the same target and conditioning information. (5) Posterior model probabilities follow Bayes’ rule from model priors and marginal likelihoods; stacking weights are selected to optimize predictive scores.