Covariance-Aware Propagation of Parameter Uncertainty

Why parameter covariance matters

A simulated exposure-response curve is only as informative as the parameter draws used to create it. If estimated parameters are correlated, drawing each parameter independently changes the combinations passed into the model. The resulting predictions can be too variable or not variable enough, and their shape can differ from predictions based on the fitted joint uncertainty.

The focus here is the covariance between a response slope and a potency parameter. We hold the exposure-response model structure fixed and propagate uncertainty in its parameter estimates. This builds on the earlier material on propagating parameter uncertainty to clinical decisions: the essential detail here is to preserve the joint distribution of the estimates when generating predictions.

Definition: For uncertain parameter estimates \(\theta_1\) and \(\theta_2\), their covariance is \(\operatorname{Cov}(\theta_1,\theta_2)=E[(\theta_1-E[\theta_1])(\theta_2-E[\theta_2])]\). Positive covariance means values above their respective means tend to occur together, as do values below their means. Negative covariance means a value above one mean tends to accompany a value below the other. The correlation is \(\rho=\operatorname{Cov}(\theta_1,\theta_2)/(\operatorname{SD}(\theta_1)\operatorname{SD}(\theta_2))\), a unitless quantity between \(-1\) and \(1\).

A slope–potency model

Consider a hypothetical exposure-response model for the response increment above a fixed baseline:

R(C;S,K)=\frac{S C}{C+K}

Here \(C\) is exposure in mg/L, \(S\) is the maximum response increment in response points, and \(K\) is the potency parameter in mg/L: at \(C=K\), the predicted increment is \(S/2\). The parameters \(S\) and \(K\) are estimated from data. For a given exposure \(C\), increasing \(S\) increases the prediction, whereas increasing \(K\) decreases it.

Suppose the estimated parameters are approximately jointly normally distributed on their stated scales, with means \(40\) response points and \(10\) mg/L, standard errors \(6\) response points and \(2\) mg/L, and correlation \(0.75\). This is an approximation to parameter-estimation uncertainty, not a claim that responses in a patient are normally distributed. The corresponding covariance matrix is:

\Sigma= \begin{pmatrix} 36 & 9\\ 9 & 4 \end{pmatrix}

The off-diagonal entry is \(0.75(6)(2)=9\), with units of response points times mg/L. It is not a correlation: using \(0.75\) in place of \(9\) in the covariance matrix would give the wrong scale. Before using such a matrix, check that it is symmetric and positive semidefinite. In this two-parameter case, its determinant is \(36(4)-9^2=63>0\), so the matrix is positive definite.

Worked Example: Generating a Correlated Parameter Draw

Generate one draw from the approximate joint distribution above. Let \(Z_1\) and \(Z_2\) be independent standard normal random variables. A construction that has the specified means, standard errors, and correlation is:

S=40+6Z_1,\qquad K=10+2\left(0.75Z_1+\sqrt{1-0.75^2}\,Z_2\right)

For the illustrative values \(Z_1=1\) and \(Z_2=0\), the draw is \(S=40+6(1)=46\) response points and \(K=10+2(0.75)(1)=11.5\) mg/L. At \(C=10\) mg/L, the predicted response increment is:

R(10;46,11.5)=\frac{46(10)}{10+11.5} =\frac{460}{21.5}=21.40\text{ response points}

The shared \(Z_1\) contributes to both parameters, which induces their positive correlation. If instead the same slope draw \(S=46\) were paired with an independent potency draw at its mean, \(K=10\) mg/L, the prediction would be \(460/20=23.00\) response points. This one-draw comparison illustrates that pairing matters; it is not, by itself, a summary of the full prediction distribution.

How covariance enters prediction uncertainty

For a differentiable prediction \(g(\boldsymbol{\theta})\), a first-order Taylor approximation around the mean parameter vector gives a local approximation to its variance. Write the gradient as \(\nabla g\) and the parameter covariance matrix as \(\Sigma\). Then:

\operatorname{Var}\{g(\boldsymbol{\theta})\} \ \approx\ \nabla g^{\mathsf T}\Sigma\nabla g

For two parameters, expanding this matrix product gives:

\operatorname{Var}\{g(\theta_1,\theta_2)\} \ \approx\ \left(\frac{\partial g}{\partial \theta_1}\right)^2 \operatorname{Var}(\theta_1) + \left(\frac{\partial g}{\partial \theta_2}\right)^2 \operatorname{Var}(\theta_2) + 2\frac{\partial g}{\partial \theta_1} \frac{\partial g}{\partial \theta_2} \operatorname{Cov}(\theta_1,\theta_2)

The last term is the covariance contribution. It follows directly from expanding \(\nabla g^{\mathsf T}\Sigma\nabla g\); it disappears only when the covariance is zero or one of the local derivatives is zero. The approximation is most reliable when parameter uncertainty is modest and the prediction is reasonably close to linear over the region of likely parameter values. It is a diagnostic and comparison tool, not a replacement for simulation when the relationship is strongly nonlinear.

Worked Example: Calculating the Covariance Contribution

At \(C=10\) mg/L and the mean estimates \(S=40\), \(K=10\), the prediction is \(R=20\) response points. Its derivatives are:

\frac{\partial R}{\partial S}=\frac{C}{C+K}=0.5, \qquad \frac{\partial R}{\partial K}=-\frac{SC}{(C+K)^2}=-1

Using the standard errors and covariance from the example above, the local variance approximation that retains covariance is:

\operatorname{Var}(R)\approx (0.5)^2(36)+(-1)^2(4)+2(0.5)(-1)(9) =9+4-9=4

The approximate standard deviation is \(\sqrt{4}=2.00\) response points. If the same marginal variances were used but the parameters were treated as independent, the covariance term would be omitted: the approximate variance would be \(9+4=13\), and the standard deviation would be \(\sqrt{13}=3.606\) response points. Here, positive covariance reduces prediction variance because the derivatives have opposite signs: a draw with a higher slope tends also to have higher potency, and those changes pull the prediction in opposing directions. This is a local approximation; exact simulated summaries need not match it precisely.

Compare joint and independent pairings

A simulation propagates parameter uncertainty by evaluating the model for many parameter draws. With joint draws, each draw keeps the slope and potency paired according to their estimated dependence. Independent draws keep the same marginal distributions but erase that pairing. The comparison below uses selected values one standard error above or below each mean. It illustrates the effect of pairing, but is deliberately not presented as a Monte Carlo estimate or as a faithful sample from the stated correlation of \(0.75\).

Worked Example: Same Marginals, Different Prediction Pairings

Use \(C=10\) mg/L and the same marginal values: \(S=46\) or \(34\) response points, and \(K=11.5\) or \(8.5\) mg/L. Pair high slope with high potency, and low with low, to illustrate the direction of positive dependence:

R(10;46,11.5)=\frac{460}{21.5}=21.40,\qquad R(10;34,8.5)=\frac{340}{18.5}=18.38

These two predictions have a range of \(21.40-18.38=3.02\) response points. Now cross-pair the same marginal values, as can happen when the parameters are independently drawn:

R(10;46,8.5)=\frac{460}{18.5}=24.86,\qquad R(10;34,11.5)=\frac{340}{21.5}=15.81

The range is now \(24.86-15.81=9.05\) response points. The marginal parameter values are unchanged, but the pairings produce substantially different predictions. In a real analysis, draw from the fitted joint parameter distribution and compare that with a deliberately independent-draw analysis using the same marginal distributions. Use enough draws to stabilize the prediction summaries. The two-point calculation explains the mechanism; it does not quantify the full prediction distribution.

Simulation checks and limitations

A reliable simulation begins with the parameterization and covariance on the same scale. For example, covariance estimated for \(S\) and \(K\) cannot be applied unchanged to \(\log(S)\) and \(\log(K)\). If the fitted model uses transformed parameters, simulate on that scale using its covariance and transform each draw back before evaluating the response model. For the raw-scale normal approximation in the examples, also check that simulated values respect the model’s parameter constraints. Normal draws can be negative; if that is not acceptably rare, use a suitable joint distribution or a transformed parameterization rather than silently discarding or truncating draws.

For each simulation, retain each joint parameter draw as a row: do not shuffle one parameter column after drawing. Check that the simulated means, standard deviations, and correlation are close to their intended values. Then inspect the resulting predictions across the exposure range of interest, including whether unusual parameter draws cause implausible curves. Monte Carlo summaries themselves have sampling error, so repeat or enlarge the simulation if key quantiles or decision-relevant predictions change noticeably with the number of draws.

A plug-in prediction at the mean estimates a central curve but does not describe parameter uncertainty. Independent parameter draws preserve marginal uncertainty but generally misrepresent joint uncertainty. A local covariance calculation is fast and helps explain direction and approximate magnitude, while joint simulation accommodates nonlinear predictions and provides prediction distributions. Its quality still depends on the fitted parameter distribution and the model assumptions. Parameter covariance alone does not account for structural uncertainty, model discrepancy, or patient-to-patient variability; those are distinct sources and should not be folded into the parameter covariance term.

Common pitfall: “The two parameters have the right standard errors, so drawing them independently preserves parameter uncertainty.” This preserves only the marginal distributions. It removes their covariance, changing the combinations passed into the response model and potentially changing prediction means, intervals, and decision-relevant probabilities. Use joint draws when the fitted estimates are correlated.

Check Your Understanding

Work through these questions, then check the solutions in the same box.

Suppose two estimates have standard errors \(4\) response points and \(1.5\) mg/L, with correlation \(-0.40\). What is their covariance, including units?
Using those standard errors, suppose the local prediction derivatives are \(0.5\) and \(-1\). Calculate the approximate prediction variance with covariance and under independence.
In the slope–potency model, why does positive covariance reduce the local prediction variance at fixed positive exposure?
What is lost if the slope and potency columns are independently shuffled after joint draws are generated?

Solutions. (1) The covariance is \((-0.40)(4)(1.5)=-2.40\) response points times mg/L. (2) With covariance, the variance is \((0.5)^2(16)+(-1)^2(2.25)+2(0.5)(-1)(-2.40)=4+2.25+2.40=8.65\) response-points squared. Under independence it is \(4+2.25=6.25\) response-points squared. (3) The slope derivative is positive and the potency derivative is negative, so their product is negative; multiplying by positive covariance makes the covariance term reduce the variance. (4) Shuffling destroys the estimated pairing, producing independent parameter pairs rather than draws from the fitted joint distribution.