Reading Skewed and Multimodal Prediction Distributions

Why the Shape of a Prediction Distribution Matters

A PBPK simulation can produce a distribution of predicted exposures by propagating variation or uncertainty in model inputs. A single average may summarize its center, but it cannot show whether values cluster symmetrically around that average, extend far into one tail, or form separate groups. Those distinctions matter when describing exposure in context: a right tail can contain unusually high predictions, while two modes may reflect different scenarios rather than one smoothly varying population.

The aim here is to interpret the distribution of a prediction, not to decide whether the PBPK model is adequate for a particular use. Begin with the quantity being predicted and its units—for example, AUC in mg·h/L or \(C_{\max}\) in mg/L. Then examine the distribution’s shape, relevant quantiles, and any scenario labels attached to the simulation draws.

Definition: Let \(Y\) denote a PBPK prediction across the simulated draws, such as \(C_{\max}\) in mg/L. Its prediction distribution describes how \(Y\) varies under the specified inputs, assumptions, scenario weights, and simulation design. The cumulative distribution function is \(F(y)=P(Y\leq y)\). The \(p\) quantile is \(q_p=\inf\{y:F(y)\geq p\}\). In a finite simulation with \(N\) equally weighted draws, the empirical CDF is \(F_N(y)=N^{-1}\sum_{i=1}^N I(Y_i\leq y)\), and the empirical quantile can be defined as the smallest observed value whose empirical CDF is at least \(p\).

Reading Density Plots and Quantiles

A density plot gives a visual summary of where simulated predictions are concentrated. A peak indicates a region of relatively high density; it does not give the probability of observing the exact value at the peak. For a continuous prediction, probabilities correspond to areas under the density over intervals. For example, the area between \(a\) and \(b\) represents \(P(a\leq Y\leq b)\).

A density estimate smooths the simulated values. Its appearance can depend on the smoothing bandwidth: too much smoothing can blur separate clusters, while too little can create small bumps that are not stable features of the distribution. Compare the plot with the raw simulation values, an empirical CDF, quantiles, and scenario labels before deciding that a peak is meaningful. The density is a diagnostic display, not proof that the prediction distribution has a particular mathematical form.

Quantiles provide numerical summaries that do not require symmetry. A central 90% prediction interval can be reported as the interval from \(q_{0.05}\) to \(q_{0.95}\); it contains the middle 90% of the distribution under the chosen convention and simulation weights. It need not be centered on the mean, and it may not include every possible prediction. Quantiles are useful when a tail is long or the distribution has distinct groups, but they do not show the shape within the interval.

Worked Example: Right-Skewed Exposure Predictions

Suppose ten hypothetical PBPK simulation draws of \(C_{\max}\), in mg/L, sorted from smallest to largest, are \(1, 1, 2, 2, 3, 4, 5, 7, 12, 23\). The long upper tail suggests right skew: most values are comparatively small, with a few much larger predictions.

The mean is the sum divided by the number of draws:

\bar{x}=\frac{1+1+2+2+3+4+5+7+12+23}{10} =\frac{60}{10}=6.00\ \text{mg/L}.

The median is the average of the fifth and sixth values, \((3+4)/2=3.50\) mg/L. The mean exceeds the median because the large upper-tail values pull the mean upward. With the empirical quantile definition above, \(q_{0.05}=1\) mg/L and \(q_{0.95}=23\) mg/L: the fifth-percentile rank is the first observation, and the 95th-percentile rank is the tenth.

To see why a symmetric summary can mislead, the sample standard deviation is \(s=\sqrt{(782-60^2/10)/9}=6.847\) mg/L, because the sum of squared values is 782. A normal-style mean plus or minus 1.96 standard deviations would be \(6.00\pm1.96(6.847)=[-7.42,19.42]\) mg/L. This is not the empirical central 90% interval, and its negative lower endpoint is not physically plausible for \(C_{\max}\). Reporting the empirical quantiles describes the simulated range more faithfully here, while the density plot would show the concentration of low values and the extended upper tail.

Multimodality and Scenario Labels

A multimodal distribution has more than one region of high density. In PBPK work, separate modes may arise when predictions combine distinct physiological or model scenarios, such as alternative formulations or different assumptions about a pathway. A mixture distribution is the weighted combination of those scenario-specific distributions. Its overall shape depends on both the distributions within scenarios and the weights assigned to them.

An aggregate density plot can suggest two peaks but cannot, by itself, explain their origin. Retaining scenario labels allows the analyst to compare conditional distributions, \(P(Y\mid S=s)\), and assess whether each mode corresponds to a particular scenario. If the simulation deliberately used equal numbers of draws per scenario for computational convenience, that allocation alone does not establish that scenarios have equal probability in the intended prediction distribution. Use the appropriate scenario weights when constructing the overall mixture.

Result: For a scenario variable \(S\) with scenario weights \(w_s=P(S=s)\), the overall prediction distribution is \(F(y)=\sum_s w_s F_s(y)\), where \(F_s(y)=P(Y\leq y\mid S=s)\), \(w_s\geq0\), and \(\sum_s w_s=1\). Its mean, when finite, is \(\mu_Y=\sum_s w_s\mu_s\), where \(\mu_s=E(Y\mid S=s)\). A mixture can have a mean in a low-density region between its modes.

Worked Example: A Mean Between Two Exposure Modes

Consider 100 hypothetical \(C_{\max}\) predictions in mg/L. Scenario A has 70 draws: 35 equal 1.5 and 35 equal 2.5. Scenario B has 30 draws: 15 equal 9 and 15 equal 11. Treat these counts as the stated mixture weights, 0.70 and 0.30, for this example.

The scenario means are \(\mu_A=(1.5+2.5)/2=2.00\) mg/L and \(\mu_B=(9+11)/2=10.0\) mg/L. Thus the mixture mean is

\mu_Y=0.70(2.00)+0.30(10.0)=4.40\ \text{mg/L}.

The low values occupy 70% of the draws and the high values occupy 30%, so the mean of 4.40 mg/L lies in the gap between the two clusters. It describes the arithmetic average under the given weights, but it does not describe a typical draw near 4.40 mg/L. The empirical median is 2.5 mg/L: positions 1–35 are 1.5, positions 36–70 are 2.5, and positions 71–100 are 9 or 11, so the 50th and 51st observations are both 2.5. A density plot with scenario-specific curves would reveal the two clusters and show which scenario contributes to each.

A symmetric interval can also span a low-density gap. The within-scenario variances are 0.25 and 1.00 \((\text{mg/L})^2\). The mixture variance is \(0.70(0.25)+0.30(1.00)+0.70(0.30)(10.0-2.00)^2=13.915\;(\text{mg/L})^2\), giving a standard deviation of 3.730 mg/L. The mean plus or minus 1.96 standard deviations is \([-2.91,11.71]\) mg/L. It includes the gap between the clusters and a negative concentration, so it is a poor description of the observed mixture.

Why the Mixture Can Be Wider Than Its Scenarios

When scenario means differ, the overall variance includes both within-scenario spread and between-scenario separation. This decomposition helps explain why a mixture can be broad even if each scenario-specific distribution is relatively concentrated.

For a finite-variance prediction \(Y\) and scenario variable \(S\), write \(\mu_s=E(Y\mid S=s)\) and \(\mu=E(Y)\). Add and subtract \(\mu_S\) inside \(Y-\mu\), then expand the square:

\begin{aligned} \operatorname{Var}(Y) &=E[(Y-\mu)^2]\\ &=E[(Y-\mu_S)^2]+E[(\mu_S-\mu)^2] +2E[(Y-\mu_S)(\mu_S-\mu)]. \end{aligned}

The cross-term is zero: conditional on \(S\), \(E(Y-\mu_S\mid S)=0\), while \(\mu_S-\mu\) is fixed. Therefore,

\operatorname{Var}(Y)=E[\operatorname{Var}(Y\mid S)] +\operatorname{Var}(E[Y\mid S]).

The first term is average within-scenario variance; the second is variance among scenario means. This identity justifies examining scenario-specific distributions alongside the mixture rather than treating one overall standard deviation as a complete description.

Worked Example: Respecting Unequal Scenario Weights

Suppose a hypothetical simulation allocates 100 draws to each of two scenarios to make their conditional distributions easy to compare. Scenario A has 50 predictions at 2 mg/L and 50 at 4 mg/L, so its mean is 3 mg/L. Scenario B has 50 at 10 mg/L and 50 at 14 mg/L, so its mean is 12 mg/L. The intended target mixture assigns weights 0.80 to A and 0.20 to B.

Pooling all 200 draws without weighting gives a mean of \((3+12)/2=7.50\) mg/L, because the computational allocation is 50:50. But the target mixture mean is

\mu_Y=0.80(3)+0.20(12)=4.80\ \text{mg/L}.

The weighted distribution assigns 40% of its probability to 2 mg/L, 40% to 4 mg/L, 10% to 10 mg/L, and 10% to 14 mg/L. Under the empirical inverse-CDF convention, its median is 4 mg/L, while its 90th percentile is 10 mg/L. These summaries reflect the specified target weights, not the equal number of draws used for computational convenience. The two scenario labels should still be shown, since the relatively small high-exposure component is easy to miss in an aggregate summary.

Choosing a Summary and Checking for Failure Modes

Different displays answer different questions. A mean and standard deviation can be useful for a roughly symmetric, unimodal distribution when the average and spread are meaningful for the decision. A median and quantile interval are often more interpretable for skewed predictions. A density plot reveals shape but depends on smoothing; an empirical CDF avoids smoothing and makes threshold probabilities and quantiles visible. Scenario-faceted density plots or separate quantile summaries help identify whether apparent modes correspond to labelled scenarios. None of these displays substitutes for stating the prediction, units, input assumptions, and scenario weights.

Check the plot against numerical summaries and the simulated draws. Look for a long one-sided tail, separated clusters, gaps, or a mean far from the median. Examine whether the apparent modes persist across reasonable density bandwidths and whether they align with scenario labels. Confirm that simulated values are finite and physically meaningful for the reported quantity, and that scenario weights used in summaries match the intended prediction. A shape that changes greatly with bandwidth or is based on very few tail draws should be described cautiously.

A common mistake is to report the mean plus or minus a fixed multiple of the standard deviation as though it were automatically a prediction interval. That reasoning assumes a shape for which the chosen limits have the intended coverage; skewness, multimodality, or physical bounds can make the resulting limits misleading or impossible. The correction is to state what the interval means, calculate quantiles from the prediction distribution when those are the desired limits, and show scenario-specific results when the mixture conceals distinct groups.

Quantiles also have limitations. A central interval can conceal multimodality, and two distributions with the same median and endpoints can have very different shapes. A density plot can reveal shape but may obscure small components. In practice, use complementary summaries: for example, a density plot with scenario labels, a median, selected quantiles, and the weights used for the overall mixture. Avoid interpreting the mean as a “typical” exposure when it lies between separated modes.

Key takeaway: Interpret a PBPK prediction distribution using its shape, quantiles, and scenario structure. A mean is one summary, not a substitute for the distribution; symmetric error bars can extend into impossible values or span a region where predictions are uncommon.

Check Your Understanding

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

What does the area under a continuous prediction density between \(a\) and \(b\) represent?
A distribution has a long upper tail and a mean above its median. What feature of the distribution can explain this difference?
Why might a mixture mean be a poor description of a typical prediction when the distribution has two separated modes?
A simulation has equal numbers of draws from two scenarios, but the target weights are 0.75 and 0.25. Should an unweighted pool be used to estimate the target mixture mean?
Name one check that helps distinguish a genuine scenario-related mode from a density-estimation artifact.

Solutions. (1) It represents \(P(a\leq Y\leq b)\), the probability that the prediction falls in that interval. (2) A small number of large predictions can pull the mean upward while leaving the median closer to the central observations. (3) The mean can fall in a low-density gap between groups, where few or no predictions occur. (4) No. Use the target weights: if the scenario means are \(\mu_A\) and \(\mu_B\), the target mean is \(0.75\mu_A+0.25\mu_B\). (5) Check whether the mode aligns with scenario labels and persists in raw values or an empirical CDF, and whether it remains under reasonable changes in density bandwidth.