Why Propagate Parameter Uncertainty by Simulation?
A fitted dose-response model gives more than a predicted response at a chosen dose: it also gives uncertainty about that prediction because the fitted parameters are uncertain. When the quantity of interest is a nonlinear function of several parameters, its uncertainty can be awkward to calculate directly. Monte Carlo propagation provides a general route: draw plausible parameter values from a specified uncertainty distribution, calculate the quantity for every draw, and summarize the resulting values.
The simulated spread is not the same as simulation error. The spread describes how much the model prediction or derived dose varies across the parameter distribution. Simulation error describes how accurately a finite set of draws estimates that spread or a summary such as its mean or quantile. Increasing the number of draws reduces simulation error; it does not reduce the underlying uncertainty in the fitted parameters.
Specify the Parameter Distribution Before Simulating
For a fitted model, a common frequentist approximation is a multivariate normal distribution centered at the estimate \(\hat{\theta}\), with covariance matrix \(\hat{\Sigma}\). This approximation represents sampling uncertainty in the estimates under the fitted model and its assumptions. It is usable only if the approximation is adequate in the relevant region and the covariance matrix is valid. A covariance matrix must be symmetric and positive semidefinite. Draws should also respect parameter constraints: for example, a parameter constrained to be positive should not be used as an unconstrained normal variable if that would generate impossible values.
If the uncertainty distribution is instead a posterior distribution, posterior draws may be propagated in the same way. The interpretation differs: posterior draws represent uncertainty conditional on the model, data, and prior, whereas a covariance-based frequentist approximation represents an estimated sampling distribution. These sources of draws are not interchangeable. In either case, record what distribution was sampled and what uncertainty it is intended to represent.
For a candidate dose \(d\), write the model’s mean response as \(g(d,\theta)\). Evaluating \(g(d,\theta^{(m)})\) for every draw gives a distribution of the mean response induced by parameter uncertainty. This is not automatically a distribution of individual patient outcomes. To represent individual outcomes, one would also need to simulate the model’s residual or outcome variability. Keeping those two layers separate prevents parameter uncertainty from being mistaken for total predictive variability.
Worked Example: Propagating Uncertainty to a Fixed-Dose Response
Suppose a hypothetical fitted model has mean response \(g(d,\theta)=\theta_0+\theta_1d\), where response is measured in units and dose \(d\) is in mg. For illustration, assume independent parameter uncertainty with \(\theta_0\sim N(2,0.5^2)\) response units and \(\theta_1\sim N(0.8,0.1^2)\) response units per mg. Consider \(d=5\) mg. These values and the simulated draws below are hypothetical.
The response at the mean parameter values is \(2+0.8(5)=6\) response units. Because this model is linear in the parameters, the exact mean response under the stated distribution is also 6. Independence gives its parameter-induced variance:
Thus the parameter-induced standard deviation is \(\sqrt{0.50}=0.7071\) response units. Now suppose four illustrative parameter draws are \((1.5,0.7)\), \((2.5,0.9)\), \((1.5,0.9)\), and \((2.5,0.7)\). Their responses at 5 mg are, respectively, \(1.5+0.7(5)=5\), \(2.5+0.9(5)=7\), \(1.5+0.9(5)=6\), and \(2.5+0.7(5)=6\) response units.
The simulated mean is \((5+7+6+6)/4=6\) response units. The sample variance is \(\{(5-6)^2+(7-6)^2+(6-6)^2+(6-6)^2\}/(4-1)=2/3\), so the simulated standard deviation is \(0.8165\) response units. The difference between \(0.8165\) and the distribution’s standard deviation \(0.7071\) is unsurprising with only four draws; the example shows how propagation works, not a suitably precise simulation.
Monte Carlo Estimates and Their Simulation Error
For independent draws from the specified parameter distribution, let \(H=h(\theta)\) have finite variance. The Monte Carlo estimate of its expected value is \(\bar{h}=M^{-1}\sum_{m=1}^{M}h^{(m)}\). Its estimated Monte Carlo standard error is \(s_h/\sqrt{M}\), where \(s_h\) is the sample standard deviation of the simulated values. This standard error measures the numerical uncertainty in estimating the expectation by a finite simulation; it is not the standard deviation of the propagated uncertainty distribution itself.
The distinction follows from the variance of an average of independent, identically distributed values. If the variance of \(H\) is \(\sigma_h^2\), then the variance of the average is:
Replacing the unknown \(\sigma_h\) by \(s_h\) gives the estimated Monte Carlo standard error. The inverse-square-root relationship is important: multiplying the number of independent draws by four approximately halves the Monte Carlo standard error. It does not change the parameter distribution being propagated.
Worked Example: Checking the Precision of a Simulated Mean
Return to the hypothetical response model at 5 mg. Its parameter-induced standard deviation is \(0.7071\) response units. If \(M=10{,}000\) independent parameter draws are propagated, the Monte Carlo standard error of the simulated mean is estimated, using that standard deviation as a planning value, by:
This says that the simulation average has a typical numerical error of about 0.0071 response units under the independent-draw assumption. It does not say that parameter uncertainty is only 0.0071; the propagated standard deviation remains about 0.7071 response units. If the target precision for the simulated mean were 0.0036 response units, the approximate required number of draws would be \((0.7071/0.0036)^2\), or about 38,580. This is a planning calculation, not a guarantee that every run will meet that error.
Propagating Parameter Draws to a Target Dose
A derived dose can be calculated for every parameter draw, provided the model defines that dose for the draw. For example, if the model’s mean response is \(g(d,\theta)\) and the target response is \(T\), solve \(g(d,\theta^{(m)})=T\) for \(d\) separately at each draw. The resulting doses approximate the distribution induced by parameter uncertainty. Some draws may yield no dose in the clinically or mathematically allowed range; those outcomes must be identified and handled according to the stated estimand, rather than silently discarded.
A nonlinear transformation can produce a skewed dose distribution even when the parameter distribution is approximately normal. Therefore, the transformed draws should be inspected and summarized in ways suited to their shape. The mean dose, standard deviation, and quantiles answer different questions; a mean and standard deviation alone may conceal asymmetry or a long tail. The earlier material on target-attaining dose distributions provides the context for interpreting such summaries.
Worked Example: Transforming Parameter Draws to Target Doses
Use the hypothetical linear model \(g(d,\theta)=\theta_0+\theta_1d\), and set the target response to \(T=6\) response units. For a positive slope \(\theta_1\), the dose attaining the target is \(D=(6-\theta_0)/\theta_1\) mg. Apply this equation to the same four illustrative draws from the fixed-dose example.
For \((\theta_0,\theta_1)=(1.5,0.7)\), \(D=(6-1.5)/0.7=6.429\) mg. For \((2.5,0.9)\), \(D=(6-2.5)/0.9=3.889\) mg. For \((1.5,0.9)\), \(D=4.5/0.9=5.000\) mg. For \((2.5,0.7)\), \(D=3.5/0.7=5.000\) mg.
The simulated mean dose is \((6.429+3.889+5.000+5.000)/4=5.079\) mg, subject to rounding. Using the unrounded doses, the sample standard deviation is approximately 1.041 mg. This four-draw summary is highly imprecise and cannot support a reliable dose interval. It does, however, show why the target-dose calculation must be repeated for every parameter draw: the transformed values vary even though all draws came from one fitted model.
How Many Draws Are Needed for a Dose Quantile?
Quantiles often matter because they describe tails of the propagated distribution. For independent draws from a continuous dose distribution with density \(f_D\) that is positive and sufficiently smooth at the true \(p\)th quantile \(q_p\), the approximate Monte Carlo standard error of the estimated quantile is:
The numerator reflects the variability in the fraction of draws below the quantile. Dividing by the density converts uncertainty in that fraction into uncertainty on the dose scale: a low density near the quantile means that a small change in cumulative probability corresponds to a larger change in dose. The approximation is unreliable when the density is near zero, the distribution has a discontinuity, or the number of effective independent draws is small.
Worked Example: Simulation Error for an Upper Dose Quantile
Suppose, purely for illustration, that the propagated target dose is approximately normal with mean 10 mg and standard deviation 2.9 mg. Its 97.5th percentile is approximately \(10+1.96(2.9)=15.684\) mg. At that percentile, the density is \(\phi(1.96)/2.9\), where the standard normal density \(\phi(1.96)\) is approximately 0.05844. Therefore, \(f_D(q_{0.975})\approx0.05844/2.9=0.02015\) per mg.
With \(M=10{,}000\) independent dose draws, the approximate Monte Carlo standard error is:
This is simulation error in estimating the 97.5th percentile, not uncertainty about the true dose parameter and not a confidence statement about a future patient. With \(M=2{,}500\), the standard error is approximately twice as large, \(0.1550\) mg; with \(M=40{,}000\), it is approximately half as large, \(0.03874\) mg. These values rely on the assumed smooth normal density and independent draws. In an actual analysis, the quantile estimate and local density should be checked using the simulated distribution rather than presumed normal.
Diagnostics, Alternatives, and Limitations
A simulation can be numerically precise while propagating an unsuitable parameter distribution or an inadequate model. Check both the simulation and the assumptions that feed it. Confirm that draws have the intended center and covariance, that the covariance matrix is valid, and that simulated parameters satisfy model constraints. Plot or tabulate the propagated responses and doses to look for skewness, extreme values, missing target solutions, and implausible predictions. Repeat the calculation with a larger number of draws or independent random-number streams; substantial changes in reported summaries indicate that the simulation is not yet stable enough.
The standard errors above assume independent draws. If draws come from a dependent simulation algorithm, such as a Markov chain, treating the raw draw count as \(M\) will generally understate Monte Carlo error. Use an appropriate effective sample size or a dependence-aware Monte Carlo standard error, and check that the chain has adequately explored the relevant distribution. Rare tail behavior can require far more draws than a central mean: a stable mean does not establish that an extreme quantile is stable.
| Approach | Useful feature | Important limitation |
|---|---|---|
| Monte Carlo propagation | Handles complicated transformations with repeated model evaluation. | Has finite simulation error and depends on the chosen parameter distribution. |
| Delta method | Can give a quick local approximation to uncertainty for a smooth transformation. | May be inaccurate for strong nonlinearity, skewness, or boundary behavior. |
| Bootstrap | Can represent sampling variability by refitting across resampled data sets. | Requires suitable resampling and repeated fitting; it answers a sampling-based question, not automatically a posterior one. |
| Profile likelihood | Can examine likelihood-based uncertainty for a derived quantity. | Requires likelihood calculations and may be more computationally involved. |
Monte Carlo propagation is especially useful when a derived quantity is nonlinear or difficult to express with a closed-form variance. It cannot repair model misspecification, weak information in the data, or an inappropriate uncertainty distribution. Nor does increasing \(M\) make a normal approximation for fitted parameters more accurate: it only estimates the consequences of that approximation more precisely. Report the distribution being propagated, the number and type of draws, the summaries, and a relevant Monte Carlo error assessment.
Check Your Understanding
Answer each question, then compare with the solutions below.
- If the propagated response values have a standard deviation of 0.6 response units, what is the approximate Monte Carlo standard error of their mean for 900 independent draws?
- Why is the standard deviation of propagated predictions not the Monte Carlo standard error of their mean?
- A fitted parameter distribution generates some negative values for a parameter that must be positive. What should be checked before interpreting the propagated results?
- If the number of independent draws is multiplied by nine, by what factor does the approximate Monte Carlo standard error change?
- Why can a dose quantile require more simulation draws than a mean?
Solutions. (1) \(0.6/\sqrt{900}=0.6/30=0.020\) response units. (2) The standard deviation describes the spread of the propagated uncertainty distribution; the standard error of the mean estimates the numerical error of the simulated average and is approximately the standard deviation divided by \(\sqrt{M}\). (3) Check whether the parameter distribution and scale are appropriate and whether the constraint should be imposed through a transformed parameter or another justified distribution; invalid draws should not be silently retained or discarded. (4) It is divided by \(\sqrt{9}=3\). (5) A quantile, particularly in a tail, depends on relatively few draws near that tail and on the local density; a low density near the quantile increases its Monte Carlo error.