Posterior Draws For Bayesian Simulation

From Posterior Draws to Predictions

A Bayesian fit represents uncertainty about model parameters with a posterior distribution. For simulation, the important practical step is to use draws from that distribution as inputs to the model, rather than reducing them immediately to one point estimate. Each draw gives a plausible parameter vector conditional on the observed data, and each vector can generate its own prediction under a specified dose, covariate profile, and outcome model.

This procedure carries parameter uncertainty into predictions. If the target is a new individual or a noisy future observation, it can also include variability generated by the model for individuals and outcomes. Those are distinct sources of variation: a posterior draw changes the model parameters, while a subject-level draw generates a new subject or outcome conditional on those parameters.

Definition: Let $\theta$ denote the full vector of model parameters and $D$ the observed data. Posterior draws $\theta^{(1)},\ldots,\theta^{(M)}$ are draws from the posterior distribution $p(\theta\mid D)$, or appropriately weighted draws representing it. For a prediction $Y$, the posterior predictive distribution is $$ p(Y\mid D)=\int p(Y\mid\theta)\,p(\theta\mid D)\,d\theta. $$ In a pharmacometric model, $p(Y\mid\theta)$ may include new-subject random effects and residual or outcome variability, depending on the prediction target.

The vector notation matters. If $\theta$ contains clearance, volume, and other parameters, the entries from one posterior draw belong together. A draw is not a collection of independent marginal values to be freely recombined. Its joint distribution carries dependence that can materially affect nonlinear quantities such as concentrations, exposure, or probability of target attainment.

Assumptions and a Draw-by-Draw Workflow

The calculation is meaningful only if the draws represent the posterior distribution for the fitted model and the conditional simulation represents the intended prediction. In particular, the model, likelihood, prior, parameterization, and data used for the fit must be those intended for the analysis. Retained draws must adequately explore the posterior, including its dependence and tails. If draws are weighted, the weights must be used in posterior summaries; treating unequally weighted draws as equally likely changes the target distribution.

The conditional simulation also needs to match the target. A prediction for a typical individual may use fixed covariates and omit new individual random effects. A prediction for a new patient may require drawing a new random effect from its population distribution. A future measured concentration may additionally require the residual-error model. These choices are part of defining the prediction, not merely technical details of the simulation.

Define the target.
Specify the quantity, time, dose or regimen, covariates, and whether the target concerns a typical individual, a new individual, or a future observed outcome.

Retain complete posterior draws.
For each retained draw, keep the full parameter vector intact. Apply any required transformations consistently with the fitted parameterization.

Run the conditional simulation.
Evaluate the prediction using that draw. Generate subject-level random effects or outcome noise only when required by the defined target.

Summarize the simulated target.
Use the collection of draw-specific predictions, with appropriate weights if applicable, to estimate posterior predictive means, quantiles, or probabilities.

The main result follows directly from integration. For equally weighted draws, replacing the integral by an average gives $$ p(Y\mid D)\ \approx\ \frac{1}{M}\sum_{m=1}^{M}p(Y\mid\theta^{(m)}), $$ where this expression describes a mixture of conditional distributions. To simulate from that mixture, select a posterior draw and then simulate $Y$ conditionally on it. A Monte Carlo estimate of a predictive mean is therefore $$ \widehat{E}(Y\mid D)=\frac{1}{M}\sum_{m=1}^{M}Y^{(m)}, $$ where $Y^{(m)}$ is simulated under $\theta^{(m)}$. If each parameter draw generates multiple subjects or outcomes, average those outcomes within and across draws according to the target being estimated.

Worked Example: Predicting Concentration From Joint Draws

Suppose three hypothetical, equally weighted posterior draws for a one-compartment concentration model are shown below. The dose is 100 mg, and the target is concentration 4 hours after a single dose, with no new-subject or residual variability added. The model is $C(t)=D/V\exp[-(CL/V)t]$, where $CL$ is in L/h and $V$ is in L.

Draw	$CL$ (L/h)	$V$ (L)
1	1	10
2	2	10
3	1	20

For draw 1, substitute its paired values into the model: $$ C^{(1)}(4)=\frac{100\ \mathrm{mg}}{10\ \mathrm{L}}\exp\left[-\frac{1\ \mathrm{L/h}}{10\ \mathrm{L}}(4\ \mathrm{h})\right] =10\exp(-0.4)=6.703\ \mathrm{mg/L}. $$

For draw 2, $CL/V=2/10=0.2\ \mathrm{h}^{-1}$, so $C^{(2)}(4)=10\exp(-0.8)=4.493\ \mathrm{mg/L}$. For draw 3, $CL/V=1/20=0.05\ \mathrm{h}^{-1}$, so $C^{(3)}(4)=5\exp(-0.2)=4.094\ \mathrm{mg/L}$. Their equally weighted average is $(6.703+4.493+4.094)/3=5.097\ \mathrm{mg/L}$.

The three pairs are the joint draws; they are not instructions to combine any clearance with any volume. The example's small number of draws is only for transparent arithmetic. An actual posterior summary needs enough representative draws for adequate Monte Carlo precision.

Why Joint Dependence and Nonlinearity Matter

The draw-by-draw calculation preserves posterior dependence automatically, provided each complete vector is kept together. Independently shuffling or resampling parameter columns changes their joint distribution, even if every column still has the same marginal values. Since pharmacometric predictions are generally nonlinear functions of parameters, the resulting prediction distribution can change.

The same issue explains why applying the model once to the posterior mean parameters is generally not equivalent to averaging predictions. In the example above, the mean clearance is $4/3$ L/h and the mean volume is $40/3$ L. The concentration calculated from those means is $$ \frac{100}{40/3}\exp\left[-\frac{(4/3)}{(40/3)}(4)\right] =7.5\exp(-0.4)=5.027\ \mathrm{mg/L}. $$ This differs from the draw-wise average of 5.097 mg/L. Neither small numerical difference should be overinterpreted here; the point is that substitution and averaging are different operations. For nonlinear targets, the difference can be consequential.

Worked Example: What Independent Column Pairing Changes

Use the same three marginal clearance values, 1, 2, and 1 L/h, and volume values, 10, 10, and 20 L. The original paired draws are $(1,10)$, $(2,10)$, and $(1,20)$. If the columns are independently rearranged into $(1,10)$, $(1,10)$, and $(2,20)$, the marginal values are unchanged, but the joint pairings have changed.

At 4 hours, the first two rearranged pairs each predict $10\exp(-0.4)=6.703\ \mathrm{mg/L}$. The third predicts $5\exp(-0.4)=3.352\ \mathrm{mg/L}$. The average is $(6.703+6.703+3.352)/3=5.586\ \mathrm{mg/L}$, rather than 5.097 mg/L from the original joint draws.

This is not a correction to the original estimate. It is a different simulation based on a different dependence structure. Retaining complete posterior rows is therefore essential when making predictions from correlated parameters.

Parameter Uncertainty Versus Subject and Outcome Variability

A posterior predictive simulation may have more than one random layer. The posterior distribution describes uncertainty about parameters conditional on the fitted data and model. Conditional on a particular parameter draw, the model can generate different subjects or outcomes through random effects and residual error. For a new individual, for example, one may draw $\eta$ from the model's random-effects distribution and then calculate that individual's outcome using $\theta^{(m)}$ and $\eta$.

The distinction can be stated with a variance decomposition. If $Y$ is a future outcome and $\theta$ is uncertain given $D$, then $$ \operatorname{Var}(Y\mid D) =E\{\operatorname{Var}(Y\mid\theta,D)\mid D\} +\operatorname{Var}\{E(Y\mid\theta,D)\mid D\}. $$ The first term is the average conditional outcome variability, including whatever subject or observation variability the model includes. The second term is variation in the conditional mean across posterior parameter draws. To justify the identity, write $Y-E(Y\mid D)=[Y-E(Y\mid\theta,D)]+[E(Y\mid\theta,D)-E(Y\mid D)]$, square, and take expectations given $D$. The cross-product has expectation zero because $E[Y-E(Y\mid\theta,D)\mid\theta,D]=0$. The two remaining squared terms give the two terms on the right.

Worked Example: Separating Parameter and Outcome Variation

Suppose a hypothetical model gives posterior probability of target attainment $p=0.40$ under one parameter draw and $p=0.60$ under another, with equal posterior weights. For each new subject, target attainment $Y$ is a Bernoulli outcome: $Y=1$ for attainment and $Y=0$ otherwise. This simplified example treats each subject outcome as conditionally independent given the draw.

The predictive mean is $(0.40+0.60)/2=0.50$. Conditional Bernoulli variance is $p(1-p)$, so the average within-draw outcome variance is $[0.40(0.60)+0.60(0.40)]/2=0.24$. The posterior variance of the conditional probability is $[(0.40-0.50)^2+(0.60-0.50)^2]/2=0.010$. Thus total predictive variance is $0.24+0.010=0.250$, which also equals $0.50(1-0.50)$.

Here 0.24 represents variation in individual binary outcomes conditional on the parameter draw, while 0.010 represents posterior variation in the target-attainment probability. A simulated data set of subjects contains both; reporting only the spread of the posterior probabilities would omit outcome variability, while reporting only variation among outcomes at one fixed probability would omit parameter uncertainty.

Diagnostics, Alternatives, and Limitations

Before interpreting a predictive distribution, check that the simulation corresponds to its stated target and that the numerical outputs support the calculation. Useful checks include:

Posterior exploration: review chain traces and relevant convergence summaries. Poor mixing or chains that explore different regions can make retained draws unrepresentative, particularly for joint dependence.
Draw integrity: verify that each row contains a complete parameter vector from one retained draw and that transformations, units, and parameter names match the simulation model.
Predictive computation: check that model execution succeeded for all intended draws and that times, doses, covariates, and output units match the target definition.
Monte Carlo precision: recognize that autocorrelated posterior draws contain less information than the same number of independent draws. Use appropriate effective-sample-size or Monte Carlo error calculations rather than assuming every retained iteration is independent.
Conditional variability: verify that random effects and residual error are included or omitted intentionally. Their absence may be right for a conditional typical-subject target and wrong for a new-subject or future-observation target.

A plug-in prediction uses one parameter estimate and is simple, but does not propagate posterior parameter uncertainty. A multivariate normal approximation or bootstrap-based parameter uncertainty propagation can provide other ways to represent parameter uncertainty; they are not automatically equivalent to posterior draws, because they arise from different inferential constructions and assumptions. Posterior simulation is especially useful when the fitted Bayesian model provides joint draws directly, but it remains conditional on the chosen model and prior. It cannot account for model misspecification merely by generating many simulations.

A common mistake is to draw each parameter independently from its posterior marginal distribution, simulate many subjects, and call the result a posterior predictive distribution. This preserves marginal distributions but usually breaks posterior dependence, so it can produce parameter combinations that do not represent the fitted posterior. The correction is to use intact joint vectors from the same retained posterior draws. Another mistake is to label the spread across posterior parameter draws as total variability among future subjects. That spread measures parameter uncertainty in the conditional prediction; total predictive variation also depends on the conditional subject and outcome model.

Key takeaway: For each posterior draw, preserve the full parameter vector, generate only the conditional variability required by the target, and summarize predictions across draws. This carries posterior uncertainty into the prediction without confusing it with variability among simulated individuals or observations.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A model has posterior draws of clearance and volume stored in paired rows. What goes wrong if the two columns are independently shuffled before concentration prediction?
Why can a prediction at the posterior mean parameter vector differ from the average of predictions across posterior draws?
A target concerns a future measured concentration in a new individual. Name two conditional sources of variability that may need to be simulated in addition to drawing model parameters.
In a predictive variance decomposition, what does the variation of conditional means across posterior draws represent?
Why should autocorrelated posterior iterations not automatically be treated as independent draws when judging Monte Carlo precision?

Solutions. (1) Shuffling changes the joint posterior dependence and thus can change predictions from nonlinear functions of clearance and volume. Keep each complete vector intact. (2) In general, applying a nonlinear function to averaged parameters is not the same operation as averaging that function's values over draws. (3) Depending on the model and target, simulate a new individual random effect and residual or observation error. (4) It represents posterior parameter uncertainty in the conditional mean prediction. (5) Correlation between iterations reduces the effective amount of independent information, so treating every iteration as independent can understate Monte Carlo error.