Why Derived Pharmacokinetic Parameters Need Their Own Uncertainty
A fitted NLME model often estimates parameters that are not the final quantities of interest. A report may need an elimination half-life, a clearance-to-volume ratio, or steady-state exposure, each calculated from fitted parameters. The uncertainty of such a quantity depends not only on the uncertainty of each fitted parameter, but also on how the function changes with those parameters and on their covariances.
The delta method provides a first-order approximation to that uncertainty. It uses the fitted parameter covariance matrix and the gradient of the derived quantity. This tutorial concerns uncertainty in a smooth function of the fitted population parameters. It does not turn that approximation into an exact description of the function’s sampling distribution.
The covariance matrix must be expressed on the same parameter scale used to calculate the gradient. For example, if a model estimates \(\log(CL)\), the gradient must be with respect to \(\log(CL)\), not \(CL\). The matrix’s off-diagonal entries matter whenever the function depends on multiple correlated parameters.
Why the Quadratic Form Gives the Approximation
The justification is a first-order Taylor expansion around the fitted value. For a parameter vector near \(\hat{\boldsymbol{\theta}}\),
If the fitted parameter vector has approximate covariance \(\widehat{\boldsymbol{\Sigma}}\), the variance of the linear term is the variance of a weighted sum. For a vector \(\mathbf{a}\), the variance of \(\mathbf{a}^{\mathsf T}\boldsymbol{\theta}\) is \(\mathbf{a}^{\mathsf T}\widehat{\boldsymbol{\Sigma}}\mathbf{a}\). Taking \(\mathbf{a}=\nabla g(\hat{\boldsymbol{\theta}})\) gives the stated formula. The approximation is therefore most credible when the function is smooth and nearly linear over the part of parameter space that contributes appreciably to estimation uncertainty.
Under the usual large-sample approximation that \(\hat{\boldsymbol{\theta}}\) is approximately multivariate normal, the linearized \(g(\hat{\boldsymbol{\theta}})\) is approximately normal as well. A two-sided 95% interval is then
This is an approximate confidence interval: under repeated sampling and the stated approximation, intervals constructed this way would cover the fixed derived parameter about 95% of the time. It does not mean there is a 95% probability that the fixed parameter lies in this particular interval.
Worked Example: Half-Life From Log Clearance and Volume
Hypothetical values. Suppose the fitted parameter vector is \((\eta,V)\), where \(\eta=\log(CL)\). At the fitted values, \(CL=5.00\) L/h and \(V=10.0\) L. The estimated covariance matrix, in the order \((\eta,V)\), is
For a one-compartment elimination rate \(k=CL/V\), the half-life is \(t_{1/2}=\log(2)/k=\log(2)V/CL\). In terms of the fitted coordinates, \(g(\eta,V)=\log(2)V\exp(-\eta)\). At the estimates, \(g(\hat{\eta},\hat{V})=\log(2)(10)/5=1.3863\) h. The gradient is
Substituting the gradient and covariance matrix into the quadratic form gives
Thus \(SE(\hat{t}_{1/2})=\sqrt{0.3075}\approx0.5545\) h. The approximate 95% interval is \(1.3863\pm1.96(0.5545)\), or approximately \([0.2995,\ 2.473]\) h. It is a local normal-approximation interval for half-life, based on the fitted covariance matrix.
The covariance contribution is negative here, reducing the variance. If the covariance were incorrectly set to zero while retaining the same marginal variances, the variance would be \(0.3075+0.3075\approx0.6150\ \mathrm{h^2}\), and the standard error would be approximately \(0.7842\) h. This difference illustrates why the full covariance matrix, rather than only its diagonal, is needed.
Other Common Pharmacokinetic Functions
For \(CL\) and \(V\), the ratio \(CL/V\) is the elimination rate constant in a one-compartment model. Its gradient with respect to \((CL,V)\) is \((1/V,\ -CL/V^2)\). For a linear pharmacokinetic model with fixed dose, steady-state AUC over one dosing interval is \(AUC_{\tau,ss}=Dose/CL\). If the fitted coordinate is \(CL\), its derivative is \(-Dose/CL^2\); if the fitted coordinate is \(\log(CL)\), its derivative is \(-Dose/CL\). These are different gradients because they refer to different parameter scales.
Worked Example: Clearance-to-Volume Ratio
Hypothetical values. Suppose \(CL=6.00\) L/h and \(V=30.0\) L, with estimated covariance matrix
where the parameter order is \((CL,V)\). The estimated ratio is \(\hat{k}=6/30=0.2000\ \mathrm{h^{-1}}\). Its gradient is
The variance estimate, including the covariance term, is
Therefore \(SE(\hat{k})\approx0.01633\ \mathrm{h^{-1}}\), and the approximate 95% interval is \(0.2000\pm1.96(0.01633)\), or \([0.1680,\ 0.2320]\ \mathrm{h^{-1}}\). For this calculation the covariance matrix is internally plausible: its determinant is \(0.36(9.00)-(1.20)^2=1.80>0\), and its diagonal entries are positive. Those checks do not establish that the fitted model is correct, but they help identify an invalid covariance input.
Worked Example: Steady-State AUC With Log Clearance
Hypothetical values. A fixed dose of 500 mg is administered at steady state under linear pharmacokinetics. Suppose the fitted clearance is \(5.00\) L/h and the estimated variance of \(\eta=\log(CL)\) is \(0.0100\). Then \(AUC_{\tau,ss}=Dose/CL=500/5=100.0\ \mathrm{mg\cdot h/L}\). As a function of \(\eta\), this is \(g(\eta)=Dose\exp(-\eta)\), with derivative \(g'(\eta)=-Dose\exp(-\eta)\). At the fitted value, the derivative is \(-100.0\ \mathrm{mg\cdot h/L}\) per unit of \(\eta\).
The delta-method variance is \((-100.0)^2(0.0100)=100.0\ (\mathrm{mg\cdot h/L})^2\), so the standard error is \(10.00\ \mathrm{mg\cdot h/L}\). The approximate 95% interval is \(100.0\pm1.96(10.00)=[80.4,\ 119.6]\ \mathrm{mg\cdot h/L}\). This calculation treats dose as fixed and propagates uncertainty in clearance only. If dose were uncertain, its variance and any covariance with clearance would also need to be included.
Assumptions, Diagnostics, and Alternatives
A delta-method interval is most defensible when the fitted parameter estimates are approximately normal on their estimation scale, the covariance matrix adequately estimates their sampling covariance, and the derived function is differentiable near the fitted values. The function should not change sharply across the uncertainty region. A narrow standard error relative to the scale over which its curvature changes supports a local approximation, but does not prove it is adequate.
Before reporting an interval, check the following:
- Parameter scale: Confirm that each derivative matches the coordinates of the covariance matrix. If a parameter was estimated on a log scale, use the chain rule on that scale.
- Covariances and ordering: Use the full matrix and verify that its parameter order matches the gradient. A matrix paired with a gradient in a different order can produce a plausible-looking but wrong result.
- Units: Check that the gradient and covariance units yield variance in the square of the derived quantity’s units, and that the standard error and interval return to the original units.
- Numerical gradient: For a complicated function, compare the analytic gradient with finite-difference derivatives at the fitted values. Large disagreement can reveal an algebra or implementation error.
- Approximation quality: Look for strong curvature, a large standard error, or a function that becomes unstable near plausible parameter values. In such cases the symmetric interval may poorly represent uncertainty.
A common mistake is to calculate a derived quantity’s standard error from the separate parameter standard errors while assuming independence. For two parameters, the variance includes \(2g_1g_2\widehat{\operatorname{Cov}}(\hat{\theta}_1,\hat{\theta}_2)\), which can either increase or decrease the result. Another mistake is to differentiate with respect to the natural parameter while inserting a covariance matrix for a transformed parameter. The correction in both cases is to write down the parameter vector, its covariance matrix, and the gradient in the same order and on the same scale before evaluating the quadratic form.
The interval also has limitations. A first-order approximation can be poor when the function is highly curved over the uncertainty region, when estimates are near a boundary, or when the fitted covariance matrix itself is unreliable. A symmetric interval can extend into values that are implausible for the quantity. Such warning signs call for scrutiny rather than an automatic reinterpretation of the interval as exact.
| Approach | What it uses | Useful consideration |
|---|---|---|
| Delta method | Gradient and fitted parameter covariance matrix | Fast and transparent for smooth functions with a credible local approximation |
| Direct Wald interval for a fitted parameter | That parameter’s estimate and standard error | Applies to the parameter itself; a derived quantity generally needs uncertainty propagation |
| Profile likelihood | Likelihood behavior as the derived quantity is constrained | Can reflect non-quadratic likelihood behavior, but requires additional fitting and calculation |
| Bootstrap | Replicate estimates and derived quantities | Can reveal skewness or irregularity, but depends on the resampling design and successful fits |
The delta method is therefore a useful first calculation and a useful diagnostic baseline, not a universal replacement for likelihood-based or resampling approaches. The uncertainty method should match the model fit’s reliability and the behavior of the particular derived quantity.
Check Your Understanding
Try each problem, then check the solutions below.
- In one sentence, state why the delta-method variance uses the full covariance matrix rather than just marginal variances.
- For \(g(CL)=Dose/CL\) with fixed dose and a covariance matrix for \(\log(CL)\), what derivative should be used?
- Suppose a function is nearly linear near the fitted estimates, but the fitted covariance matrix has a negative diagonal entry. Is the delta-method interval ready to report? Explain.
- For a one-compartment model with \(CL=8\) L/h and \(V=40\) L, calculate \(CL/V\) and give its units.
Solutions. 1. The function’s uncertainty depends on how the parameter estimates vary together, represented by the off-diagonal covariances in the quadratic form. 2. With \(\eta=\log(CL)\), \(g(\eta)=Dose\exp(-\eta)\), so \(g'(\eta)=-Dose\exp(-\eta)=-g(\eta)\). 3. No. A covariance matrix must have nonnegative variances on its diagonal; a negative entry signals an invalid or incorrectly supplied matrix that must be investigated before computing the interval. 4. \(CL/V=8/40=0.2000\ \mathrm{h^{-1}}\), since \((\mathrm{L/h})/\mathrm{L}=\mathrm{h^{-1}}\).