Why Parameter Covariance Must Travel With the Estimates
A model handoff often transfers more than a point prediction. A downstream model may need fitted parameter estimates from an upstream model, such as clearance and volume from a population PK analysis, or an exposure-response intercept and slope. Those estimates are uncertain, and their uncertainties are often related. If one parameter tends to be estimated high when another is estimated low, or if both tend to move together, that dependence affects functions of the parameters used by the receiving model.
The central practical requirement is therefore to transfer the joint parameter uncertainty, not just each parameter's standard error. The covariance matrix records both the marginal variances and the pairwise covariances. Discarding its off-diagonal entries changes the uncertainty assigned to many downstream predictions; the direction and size of the change depend on how the prediction responds to each parameter.
Covariance is not correlation. Covariance depends on the units of both parameters: \(\operatorname{Cov}(\hat{\theta}_i,\hat{\theta}_j)\) has units equal to the product of their units. Correlation is the dimensionless covariance divided by the product of their standard errors. A handoff should retain covariance in the parameter units and parameterization in which it is used, or explicitly transform it when those change.
A covariance matrix is symmetric and positive semidefinite. In particular, a two-parameter covariance \(c\) must satisfy \(c^2 \leq \operatorname{Var}(\hat{\theta}_1)\operatorname{Var}(\hat{\theta}_2)\). A matrix that violates this condition cannot be a valid covariance matrix and signals a transcription, ordering, rounding, or estimation problem.
Propagating Covariance to a Downstream Quantity
Suppose the downstream quantity of interest is a smooth function \(g(\boldsymbol{\theta})\) of the transferred parameters. Near the fitted values, a first-order Taylor approximation gives the change in the quantity as a weighted sum of parameter changes. The weights are the derivatives of the function at the fitted values.
Taking the variance of that linear approximation gives the delta-method variance. This is the first-order uncertainty calculation used in the examples below.
For two parameters, write the gradient as \((g_1,g_2)^\mathsf{T}\), where \(g_i\) is the derivative with respect to \(\theta_i\). Expanding the matrix product shows exactly where covariance enters:
The covariance contribution can increase or decrease the result. If the two derivatives have the same sign, positive covariance increases the propagated variance, while negative covariance decreases it. If the derivatives have opposite signs, those directions reverse. A diagonal-only calculation removes the final term; it is not automatically conservative.
Worked Example: A PK Prediction With Negative Parameter Covariance
Suppose a fitted exposure-response relationship predicts mean response \(m(x)=\beta_0+\beta_1x\), with a dimensionless exposure variable \(x\). In hypothetical data, the fitted intercept is \(10\) response units and the fitted slope is \(2\) response units per unit of \(x\). The estimated parameter covariance matrix, with the intercept first and slope second, is:
The matrix is valid: its determinant is \(4(0.25)-(-0.8)^2=0.36>0\), and its diagonal entries are positive. At \(x=3\), the predicted mean response is \(10+2(3)=16\) response units. Its gradient with respect to \((\beta_0,\beta_1)\) is \((1,3)^\mathsf{T}\), so the variance due to parameter uncertainty is:
Thus the standard error of the estimated mean response is \(\sqrt{1.45}=1.204\) response units. Ignoring covariance would instead give variance \(4+2.25=6.25\), and standard error \(2.500\) response units. Under a normal approximation, the corresponding approximate 95% interval for the mean response is \(16\pm1.96(1.204)\), or \(13.64\) to \(18.36\) response units. The diagonal-only calculation gives \(16\pm1.96(2.500)\), or \(11.10\) to \(20.90\). Here, dropping the negative covariance substantially widens the interval.
These are intervals for the model's mean response, reflecting parameter uncertainty. If the target is a future individual response and the residual variance is known to be \(9\) response-units-squared and independent of parameter estimation, it must also be included. The full-covariance predictive variance is then \(1.45+9=10.45\), rather than \(1.45\); the diagonal-only version is \(6.25+9=15.25\). The two targets—mean response and a future observation—must not be confused.
Scale Changes and Linked Model Components
Model parameters are frequently reported on transformed scales. For example, a model may be fitted using log clearance and log volume, while the receiving component expects clearance and volume in natural units. Transforming point estimates alone is insufficient: the covariance matrix must be transformed as well.
For a differentiable transformation \(\boldsymbol{\phi}=q(\boldsymbol{\theta})\), first-order covariance transformation uses the Jacobian matrix \(J\), whose entries are the derivatives of each transformed parameter with respect to each original parameter:
This expression follows by applying the same first-order variance argument to the vector transformation: locally, changes in \(\boldsymbol{\phi}\) are approximately \(J\) times changes in \(\boldsymbol{\theta}\). For exponentiation of two log parameters, the Jacobian is diagonal, with the corresponding natural-scale parameter values on its diagonal. The off-diagonal covariance is multiplied by both scale factors.
Worked Example: Transforming Log-Scale PK Covariance
Suppose a hypothetical PK fit reports \(\eta_{CL}=\log(CL)=\log(5)\) and \(\eta_V=\log(V)=\log(50)\), where \(CL\) is in L/h and \(V\) is in L. Their estimated covariance matrix is:
The receiving component uses natural-scale \(CL\) and \(V\). The Jacobian of \((CL,V)=(\exp(\eta_{CL}),\exp(\eta_V))\) at the estimates is \(\operatorname{diag}(5,50)\). Therefore:
The approximate natural-scale variances are \(1\ (\mathrm{L/h})^2\) for clearance and \(225\ \mathrm{L}^2\) for volume; their covariance is \(7.5\ \mathrm{L^2/h}\). Now consider the rate \(k=CL/V\), in \(\mathrm{h}^{-1}\), whose estimate is \(5/50=0.1\ \mathrm{h}^{-1}\). On the log scale, \(\log(k)=\eta_{CL}-\eta_V\), so its gradient is \((1,-1)^\mathsf{T}\). The approximate variance of \(\log(k)\) is \(0.04+0.09-2(0.03)=0.07\). Consequently, the approximate variance of \(k\) on the natural scale is \(0.1^2(0.07)=0.0007\ (\mathrm{h}^{-1})^2\), with standard error \(0.02646\ \mathrm{h}^{-1}\).
If the log-scale covariance were discarded, the corresponding log-scale variance would be \(0.04+0.09=0.13\), giving an approximate natural-scale variance of \(0.0013\) and standard error \(0.03606\ \mathrm{h}^{-1}\). The positive covariance reduces uncertainty in the difference of log parameters. This calculation also illustrates why scale, units, and parameter definition must accompany the matrix.
When a downstream quantity depends on estimates from more than one fitted component, use a joint parameter vector. For example, let \(\boldsymbol{u}=(\boldsymbol{\theta}_A^\mathsf{T},\boldsymbol{\theta}_B^\mathsf{T})^\mathsf{T}\). Its covariance has a block form: the diagonal blocks describe uncertainty within each component, and the off-diagonal blocks describe covariance between component estimates. For a downstream function \(h(\boldsymbol{u})\), the same gradient formula applies to the full matrix. Setting cross-component covariance to zero is justified only when the estimation structure and shared information support that assumption; separate files or separate fits do not establish independence.
Worked Example: Positive Covariance Increases Response Uncertainty
A receiving exposure-response model uses a fitted intercept \(\alpha=1.0\) response unit and slope \(\beta=0.40\) response units per exposure unit. At exposure \(x=2\), its predicted mean is \(\alpha+\beta x=1.80\) response units. Suppose the estimated variances are \(0.09\) and \(0.01\), respectively, and their covariance is \(0.025\). The covariance matrix is valid because its determinant is \(0.09(0.01)-0.025^2=0.000275>0\).
The prediction gradient is \((1,2)^\mathsf{T}\). Its parameter-uncertainty variance, retaining covariance, is:
The standard error is \(\sqrt{0.23}=0.4796\) response units, giving an approximate 95% interval for the mean response of \(1.80\pm1.96(0.4796)\), or \(0.8600\) to \(2.740\) response units. A diagonal-only calculation gives variance \(0.13\), standard error \(0.3606\), and interval \(1.093\) to \(2.507\). Since both derivatives are positive and the covariance is positive, removing the covariance contribution understates uncertainty.
Practical Handoff Checks and Failure Modes
A sound handoff makes the covariance usable and auditable. The receiving analyst needs more than a matrix of numbers: parameter names and ordering, units, scale or transformation, the fitted values to which the covariance corresponds, and the meaning of the reported uncertainty. If a receiving model changes parameterization, transform the covariance before using it. If a parameter is fixed, re-estimated, or defined differently in the receiving model, document that change rather than silently treating the estimates as identical.
- Check ordering. A matrix with the right numbers in the wrong parameter order gives incorrect propagated uncertainty. Match every row and column to the parameter vector.
- Check validity and units. Confirm symmetry, positive semidefiniteness, and units. Small negative eigenvalues may result from rounding, but a materially invalid matrix requires investigation.
- Check the target. State whether the calculation concerns a model mean, an individual future outcome, or another derived quantity. Add residual variability only when the target requires it and the model assumptions justify the addition.
- Check sensitivity. Compare full-covariance and diagonal-only calculations as a diagnostic, not as two interchangeable methods. Large differences signal that dependence matters for this prediction.
- Check approximation quality. The delta method is local. If the function is strongly nonlinear over the plausible parameter range, its first-order interval can be misleading even when the covariance was transferred correctly.
A common mistake is to assume that omitting covariance is conservative because it “leaves out” an uncertainty contribution. The examples show why that reasoning fails: the contribution is \(2g_1g_2\operatorname{Cov}(\hat{\theta}_1,\hat{\theta}_2)\), which can be positive or negative. The correction is to use the full matrix and evaluate the covariance term with the actual prediction gradients.
The first-order calculation is an efficient alternative when a smooth prediction and a reliable covariance estimate are available. A joint refit may be preferable when component parameters are estimated together or when the modular handoff obscures shared-data dependence. Joint parameter draws or bootstrap refits can retain nonlinear effects, but they must preserve the joint dependence and the model structure; drawing each parameter independently from its marginal distribution does not do so. These alternatives do not remove the need to define the handoff parameters, scales, and target quantity precisely.
Check Your Understanding
Answer each question, then compare your reasoning with the solution provided.
- A prediction is \(g=\theta_1+\theta_2\). The two parameter variances are \(1\) and \(4\), and their covariance is \(1.5\). What is the approximate prediction variance, and what would a diagonal-only calculation give?
Solution: The gradient is \((1,1)^\mathsf{T}\). The full variance is \(1+4+2(1.5)=8\). The diagonal-only variance is \(1+4=5\); it understates uncertainty because both derivatives and the covariance are positive. - For \(g=\theta_1-\theta_2\), the parameter covariance is positive. Does the covariance term increase or decrease the propagated variance?
Solution: The gradient components have opposite signs, so \(2g_1g_2\operatorname{Cov}(\hat{\theta}_1,\hat{\theta}_2)\) is negative. Positive covariance decreases the first-order propagated variance for this difference. - Why is it not enough to exponentiate fitted log-scale parameters before sending them to a receiving model?
Solution: The covariance must also be transformed. For a differentiable transformation, use \(\widehat{\Sigma}_{\phi}\approx J\widehat{\Sigma}_{\theta}J^\mathsf{T}\); otherwise the receiving model gets uncertainty on the wrong scale. - A model predicts a mean response, but the analyst reports an interval for a future individual response. What additional uncertainty may be needed?
Solution: If the model includes residual outcome variability, and the assumptions support its addition independently of parameter estimation, include that residual variance as well as the parameter-uncertainty variance. State clearly that the target is a future observation rather than the mean response.