Propagating PK Uncertainty With the Delta Method

Why PK uncertainty needs to reach the PD estimate

In a sequential analysis, a PD model may be fitted using concentrations predicted from estimated PK parameters. Treating those predictions as known can make the PD estimate appear more precise than it is: the predictions depend on PK parameters that were themselves estimated with uncertainty. The delta method provides a first-order approximation to how that PK uncertainty changes a PD estimate.

The calculation is local. It uses the slope of the PD estimate as a function of the PK parameters, together with the covariance among those PK parameter estimates. A PD estimate that is highly sensitive to an uncertain PK parameter can inherit substantial uncertainty; a parameter with little sensitivity may inherit little. Correlations among PK estimates matter because their effects can reinforce or offset each other.

Definition: Let $\hat{\eta}$ be a vector of estimated PK parameters, with estimated covariance matrix $V_{\eta}$. Let $\hat{\beta}=h(\hat{\eta})$ be a scalar PD estimate that depends on those PK estimates. The delta-method gradient is $g=\nabla_{\eta}h(\hat{\eta})$, the vector of partial derivatives of $h$ with respect to the PK parameters. The first-order approximation to the variance contributed by PK uncertainty is $$ \operatorname{Var}_{\mathrm{PK}}(\hat{\beta})\approx g^{T}V_{\eta}g. $$ For a vector of PD estimates $\hat{\boldsymbol{\beta}}=h(\hat{\eta})$, let $J$ be the matrix of partial derivatives, with one row per PD estimate and one column per PK parameter. Then the corresponding covariance contribution is approximately $JV_{\eta}J^{T}$.

This is a propagation term, not necessarily the entire uncertainty of the PD estimate. Suppose the PD fitting procedure also has conditional variance $V_{\beta\mid\eta}$ when the PK parameters are held fixed. If the PK and conditional PD estimation errors are independent under the staged analysis, a first-order total covariance approximation is

V_{\beta}\approx V_{\beta\mid\eta}+JV_{\eta}J^{T}.

For one PD estimate, the standard error is the square root of the corresponding total variance. The independence assumption is important: if PK and PD estimation errors are correlated, the variance can also contain cross-covariance terms. The simple sum is not justified merely because the analysis was performed in stages.

Derivation and interpretation

The approximation follows from a first-order Taylor expansion around the fitted PK parameter vector. For a scalar PD estimate, write the estimation error in PK parameters as $\hat{\eta}-\eta$. Near $\eta$,

h(\hat{\eta})\approx h(\eta)+g^{T}(\hat{\eta}-\eta).

Taking the variance of the linear term gives

\operatorname{Var}\{g^{T}(\hat{\eta}-\eta)\} =g^{T}\operatorname{Var}(\hat{\eta})g \approx g^{T}V_{\eta}g.

For multiple PD estimates, the same expansion is $\hat{\boldsymbol{\beta}}\approx h(\eta)+J(\hat{\eta}-\eta)$. Taking its covariance gives $JV_{\eta}J^{T}$. This derivation also explains the method's limitation: it replaces a possibly curved function by its tangent near the estimate. The approximation is most credible when PK uncertainty is not so large that the function changes curvature substantially over plausible parameter values.

The diagonal entries of $V_{\eta}$ alone are not enough when PK estimates are correlated. In the scalar expression, each variance contribution is weighted by a squared derivative, while each covariance contributes a cross-term. Depending on the signs of the derivatives and covariance, that cross-term can increase or decrease the propagated variance. The full covariance matrix must therefore be used.

Worked calculations

Worked Example: A Linear PD Slope From an Uncertain Concentration

Consider a hypothetical one-compartment prediction $C(t)=A\exp(-kt)$, with fitted PK parameters $\eta=(\log A,k)$. At $t=2$ h, suppose $\hat{A}=10.00$ mg/L and $\hat{k}=0.2000$ h$^{-1}$. A PD response contrast of $13.406$ response units is divided by the predicted concentration to estimate a linear slope $\hat{\beta}=\Delta E/C(t)$. The fitted PK covariance matrix is $$ V_{\eta}= \begin{pmatrix} 0.0100 & -0.001000\\ -0.001000 & 0.0004000 \end{pmatrix}. $$

First calculate the concentration and slope: $$ \hat{C}(2)=10.00\exp(-0.2000\times2)=10.00\exp(-0.4000)=6.703\ \text{mg/L}, \qquad \hat{\beta}=\frac{13.406}{6.703}=2.000\ \frac{\text{response units}}{\text{mg/L}}. $$

Since $\hat{\beta}=\Delta E\exp(kt)/A$, its gradient with respect to $(\log A,k)$ is $(-\hat{\beta},t\hat{\beta})=(-2.000,4.000)$. The PK contribution to its variance is $$ g^{T}V_{\eta}g =(-2.000)^2(0.0100)+(4.000)^2(0.0004000) +2(-2.000)(4.000)(-0.001000) =0.06240. $$ The PK-propagated standard error is $\sqrt{0.06240}=0.2498$ response units per mg/L.

Suppose the PD estimator has conditional standard error $0.3000$ response units per mg/L when PK parameters are fixed, and assume that conditional estimation error is independent of PK estimation error. Then $$ SE(\hat{\beta})\approx\sqrt{0.06240+0.3000^2} =\sqrt{0.1524}=0.3904\ \frac{\text{response units}}{\text{mg/L}}. $$ A normal-approximation 95% interval is $2.000\pm1.96(0.3904)$, or approximately $1.235$ to $2.765$ response units per mg/L. This interval uses both stated variance contributions; it is not the conditional interval obtained by treating the PK fit as exact.

The positive cross-term arises because the gradient components have opposite signs and the PK parameter covariance is negative. Omitting that covariance would give a propagated variance of $0.04640$, rather than $0.06240$.

Worked Example: Propagating Uncertainty to an EC50 Estimate

Suppose, hypothetically, that a response $E=40.00$ units is observed at $t=2$ h, with known $E_{\max}=100.0$ units. Under the Emax relationship $E=E_{\max}C/(EC_{50}+C)$, solving for potency gives $\widehat{EC}_{50}=C(E_{\max}-E)/E$. Use the same fitted concentration and PK covariance matrix as in the preceding calculation. Assume the conditional standard error of this PD estimate, given PK, is $0.8000$ mg/L and is independent of PK estimation error.

The multiplier is $(100.0-40.00)/40.00=1.500$, so $$ \widehat{EC}_{50}=6.703(1.500)=10.05\ \text{mg/L}. $$ With $q=(E_{\max}-E)/E=1.500$, the estimate is $qA\exp(-kt)$. Its gradient with respect to $(\log A,k)$ is $(\widehat{EC}_{50},-t\widehat{EC}_{50})=(10.05,-20.11)$ in the corresponding parameter units.

Using unrounded values, the propagated variance is $$ (10.0548)^2(0.0100)+(-20.1096)^2(0.0004000) +2(10.0548)(-20.1096)(-0.001000) \approx1.577\ (\text{mg/L})^2. $$ Thus the PK-propagated standard error is $\sqrt{1.577}=1.256$ mg/L. Adding the conditional variance gives $$ SE(\widehat{EC}_{50})\approx\sqrt{1.577+0.8000^2} =1.489\ \text{mg/L}. $$ The sign pattern again makes the covariance cross-term positive. This calculation treats $E_{\max}$ and the response contrast as fixed for the propagation step; if they are estimated with relevant uncertainty, their uncertainty must also be represented in the calculation.

Worked Example: A PD Slope Based on AUC

For a hypothetical intravenous dose $D=100.0$ mg in a one-compartment model, let $AUC=D/CL$, where $CL$ is clearance. A PD slope is estimated as $\hat{\gamma}=\Delta E/AUC$, with a response contrast $\Delta E=20.00$ units. Suppose $\widehat{CL}=10.00$ L/h and $\operatorname{Var}(\widehat{\log CL})=0.09000$. The conditional PD standard error is $0.5000$ response units per mg·h/L, independent of the PK error.

The estimated AUC and slope are $$ \widehat{AUC}=\frac{100.0\ \text{mg}}{10.00\ \text{L/h}} =10.00\ \text{mg·h/L}, \qquad \hat{\gamma}=\frac{20.00}{10.00} =2.000\ \frac{\text{response units}}{\text{mg·h/L}}. $$ Since $\hat{\gamma}=\Delta E\,CL/D$, its derivative with respect to $\log CL$ is $\hat{\gamma}$. The propagated variance and standard error are therefore $$ \operatorname{Var}_{\mathrm{PK}}(\hat{\gamma}) \approx(2.000)^2(0.09000)=0.3600, \qquad SE_{\mathrm{PK}}(\hat{\gamma})=\sqrt{0.3600}=0.6000 \ \frac{\text{response units}}{\text{mg·h/L}}. $$ Combining the stated conditional variance gives a total standard error of $\sqrt{0.3600+0.5000^2}=0.7810$ response units per mg·h/L.

Assumptions, diagnostics, and alternatives

The calculation requires a differentiable mapping from PK estimates to the PD estimate in the region of interest, a suitable estimate of the PK covariance matrix, and a gradient evaluated at the fitted PK parameters. The usual total-variance formula additionally requires the specified conditional PD variance and the required independence (or a fuller covariance calculation when errors are dependent). The examples used fitted covariance matrices as given; in practice, the covariance estimate must correspond to the PK parameters in the same parameterization used for the gradient.

Check the computation before interpreting its result. Confirm parameter ordering in both the gradient and covariance matrix; a mismatch can produce plausible-looking but incorrect numbers. Confirm units for every derivative and resulting variance. Check that the covariance matrix is symmetric and positive semidefinite up to numerical tolerance. When derivatives are obtained numerically, compare them with analytic derivatives where possible, or check that small changes in the finite-difference step do not materially change the result. A negative calculated variance usually indicates an implementation error or an invalid covariance estimate, not a negative amount of uncertainty.

Large gradients, high PK uncertainty, strong parameter correlations, and a highly curved mapping are warnings that the first-order approximation may be poor. A useful diagnostic is to examine the mapping over a plausible range of PK parameter values: if small changes around the estimate produce strongly asymmetric changes in the PD estimate, a symmetric standard error may summarize uncertainty badly. Weak identification or estimates near parameter boundaries can also make a local normal approximation unreliable.

The delta method is computationally efficient and directly shows which PK parameters contribute to PD uncertainty. A conditional PD standard error alone is simpler but answers a narrower question: it treats the PK fit as fixed. A joint likelihood can represent PK and PD data together, while a resampling or simulation approach can evaluate the transformed uncertainty without relying on a first-order Taylor approximation. These alternatives have their own assumptions and computational costs; the delta method is most useful when its local approximation is defensible and the covariance inputs are available.

Common pitfall: A small conditional PD standard error does not establish that the PD estimate is precise overall. It may only reflect PD uncertainty conditional on fixed PK predictions. Include the gradient-weighted PK covariance contribution, and include cross-covariance terms when the estimation errors are not independent.

Check the calculation before relying on it

Use the questions below to test whether you can identify the relevant gradient, covariance contribution, and assumptions. Solutions are included so each calculation can be checked directly.

Check Your Understanding

Try each question first, then compare your reasoning with the solutions.

A scalar PD estimate has gradient $g=(2,1)$, and the PK covariance matrix is $\begin{pmatrix}0.04&-0.01\\-0.01&0.09\end{pmatrix}$. What is its PK-propagated variance?
If a PD estimate is $h(\hat{\eta})=5/\hat{C}$ and $\hat{C}=4.00$ mg/L, what is its derivative with respect to $C$, and what does the sign mean?
Why can two PK parameter variances be insufficient to calculate the propagated variance of a PD estimate?
A conditional PD variance is $0.16$, and the PK-propagated variance is $0.09$. Under the independence assumption, what is the total standard error?

Solutions. (1) $g^{T}V_{\eta}g=2^2(0.04)+1^2(0.09)+2(2)(1)(-0.01)=0.21$. (2) $dh/dC=-5/C^2=-5/16=-0.3125$ per mg/L at $C=4.00$ mg/L; increasing concentration decreases this estimate. (3) The covariance between PK parameter estimates contributes a cross-term, which may increase or decrease the result. (4) The total variance is $0.16+0.09=0.25$, so the total standard error is $\sqrt{0.25}=0.500$ in the units of the PD estimate.