Tutorials › Optimal Design and Model-Informed Dosing › Residual Error Models and Optimal Sampling Times

Optimal Sampling for Individual PK Models · Tutorial 51 of 497

Residual Error Models and Optimal Sampling Times

Learn to incorporate each residual-error model into a local information matrix and compare sampling schedules without treating all concentration measurements as equally precise.

Advanced 13 min read

What You'll Learn

  • Distinguish additive, proportional, and combined residual variance in concentration data
  • Derive the normal-model Fisher information when both mean and variance depend on PK parameters
  • Calculate local D-information for candidate sampling schedules
  • Understand when variance-derivative information matters and when weighted sensitivities alone are incomplete
  • Diagnose schedules that appear informative only because an error model or nominal value was misspecified
  • Compare schedules under a common parameterization and fixed residual-error assumptions

Why the Residual-Error Model Changes the Design

A concentration sample is informative about pharmacokinetic parameters partly because the predicted concentration changes with those parameters, and partly because the observation is measured with finite error. If two samples have similar sensitivities but very different residual variances, they do not contribute equally to parameter precision. Additive error gives roughly constant variance across concentrations; proportional error gives variance that grows with the predicted concentration; combined error includes both components. These distinctions can alter which sampling times are most useful.

The earlier material on D-optimal sampling establishes the use of the determinant of the information matrix as a local measure of joint parameter information. Here the key refinement is to calculate information under an explicit residual-error distribution. The design is local: predictions and their sensitivities are evaluated at specified nominal PK parameters. Residual-error components are treated as known for the comparison, and candidate schedules have the same number of samples.

Definition: For a predicted concentration \(f(t;\theta)\), additive residual error has variance \(v(t;\theta)=\sigma_a^2\); proportional residual error has variance \(v(t;\theta)=\sigma_p^2 f(t;\theta)^2\); and combined residual error has variance \(v(t;\theta)=\sigma_a^2+\sigma_p^2 f(t;\theta)^2\). Here \(\sigma_a\) is the additive standard deviation in concentration units and \(\sigma_p\) is a dimensionless proportional-error coefficient. In each case, \(Y(t)\mid\theta\) is assumed normally distributed with mean \(f(t;\theta)\) and the stated variance.

The model assumes conditionally independent observations at the scheduled times, a differentiable mean and variance, and positive variance at each time used. For a local design comparison, the nominal parameter values and residual-error coefficients are held fixed. The normal likelihood formula below also assumes the stated Gaussian residual distribution is appropriate; robust or non-Gaussian errors require a different likelihood-based information calculation.

Information When Variance Depends on the PK Parameters

Write \(f_j=\partial f/\partial\theta_j\) and \(v_j=\partial v/\partial\theta_j\). For one normal observation with mean \(f\) and variance \(v\), its Fisher information matrix has entries

$$ I_{jk}(t;\theta)=\frac{f_j f_k}{v}+\frac{1}{2}\frac{v_jv_k}{v^2}. $$

The first term is information carried by changes in the mean, reduced according to the observation variance. The second term is information carried by changes in the variance itself. It is zero for additive error when \(\sigma_a\) is fixed, but generally nonzero for proportional and combined error because their variances depend on the predicted concentration.

A Short Derivation and a Useful Special Case

For one observation \(y\), the log likelihood, omitting no parameter-dependent terms, is

$$ \ell(\theta)=-\frac{1}{2}\log(2\pi v)-\frac{(y-f)^2}{2v}. $$

Let \(\varepsilon=y-f\). Differentiating gives the score component

$$ \frac{\partial\ell}{\partial\theta_j} =\frac{f_j\varepsilon}{v} +\frac{v_j}{2v^2}(\varepsilon^2-v). $$

Under the model, \(E(\varepsilon)=0\), \(E(\varepsilon^2)=v\), \(E(\varepsilon^3)=0\), and \(E\{(\varepsilon^2-v)^2\}=2v^2\). Thus the cross-product between the two score terms has expectation zero. Taking the expected product of the scores leaves \(f_jf_k/v\) from the mean term and \(v_jv_k/(2v^2)\) from the variance term, yielding the information formula above. Independent observations contribute additively, so the schedule information is \(I(\theta)=\sum_i I(t_i;\theta)\).

A particularly convenient case is a model with log-amplitude and rate parameters, \(\theta=(\log A,k)\), and mean \(f(t)=A e^{-kt}\). Its mean sensitivity vector is \(f(t)(1,-t)^\mathsf{T}\). Under combined error, define \(r(t)=\sigma_p^2 f(t)^2/v(t)\). Since the variance derivatives are \(v_{\log A}=2\sigma_p^2 f^2\) and \(v_k=-2t\sigma_p^2 f^2\), the one-observation information is

$$ I(t;\theta)=w(t) \begin{pmatrix} 1 & -t\\ -t & t^2 \end{pmatrix}, \qquad w(t)=\frac{f(t)^2}{v(t)}+2r(t)^2. $$

The second part of \(w(t)\) is the variance information. For additive error, \(w(t)=f(t)^2/\sigma_a^2\). For proportional error with positive \(f(t)\), \(f(t)^2/v(t)=1/\sigma_p^2\) and \(r(t)=1\), so \(w(t)=1/\sigma_p^2+2\). This special structure makes pairwise schedule comparisons transparent.

Worked Example: Deriving the Weight at a Candidate Time

Consider the hypothetical one-compartment mean \(f(t)=10e^{-t/2}\), with time in hours and concentration in mg/L. Suppose the combined-error coefficients are \(\sigma_a=0.5\) mg/L and \(\sigma_p=0.1\). At \(t=2\) hours, \(f=10e^{-1}=3.6788\) mg/L and \(f^2=100e^{-2}=13.5335\) (mg/L)\(^2\).

The combined variance is \(v=0.25+0.01(13.5335)=0.385335\) (mg/L)\(^2\). Therefore \(r=0.01(13.5335)/0.385335=0.3513\). The mean contribution is \(f^2/v=13.5335/0.385335=35.12\), and the variance contribution is \(2r^2=2(0.3513)^2=0.2468\). The resulting local weight is \(w(2)=35.37\), to four significant figures.

Ignoring the variance derivative would retain only 35.12. That may be a useful mean-sensitivity approximation, but it is not the full Fisher information for the stated Gaussian combined-error model.

Comparing Schedules With D-Information

For \(n\) observations, the local D-criterion is \(\det I(\theta)\); a larger determinant indicates greater joint information about the parameter vector, for the same parameterization and parameter dimension. In the two-parameter exponential example, two samples at \(t_i\) and \(t_j\) give

$$ \det I=w(t_i)w(t_j)(t_j-t_i)^2. $$

To see this, the information matrix is the sum of \(w(t_i)(1,-t_i)^\mathsf{T}(1,-t_i)\) and \(w(t_j)(1,-t_j)^\mathsf{T}(1,-t_j)\). Expanding its two-by-two determinant cancels the terms involving \(t_i^2\) and \(t_j^2\) separately, leaving \(w(t_i)w(t_j)(t_j-t_i)^2\). Thus the schedule benefits from both large individual weights and separation in time. Separation alone is not enough: the later sample may have little weight under additive or combined error.

Worked Example: Error Model Can Change the Preferred Pair

Use the hypothetical model \(f(t)=10e^{-t/2}\) with \(\theta=(\log A,k)\), now evaluated at the nominal values \(A=10\) mg/L and \(k=1\) h\(^{-1}\). Suppose the candidate times are 0, 1, and 4 hours, and the schedule must use exactly two of them. Use \(\sigma_a=0.5\) mg/L and \(\sigma_p=0.1\). The calculated weights are as follows:

Error model\(w(0)\)\(w(1)\)\(w(4)\)
Additive400.0147.27.326
Proportional102.0102.0102.0
Combined81.2860.256.835

For additive error, the weights are \(f^2/0.25\), giving \(400.0\), \(400e^{-1}=147.1518\), and \(400e^{-4}=7.3263\). For proportional error, each weight is \(1/0.01+2=102\). For combined error at \(t=1\), \(f^2=100e^{-1}\), \(v=0.25+0.01f^2\), and \(r=0.01f^2/v\approx0.595390\). Hence \(w(1)\approx59.5390+2(0.595390)^2=60.2480\). At \(t=4\), \(r\approx0.068262\), so \(w(4)\approx6.826154+0.009319=6.83547\). At \(t=0\), \(w(0)=80+2(0.8)^2=81.28\).

For the additive model, the determinants for pairs \((0,1)\), \((0,4)\), and \((1,4)\) are respectively \(400(147.1518)(1)^2=58,861\), \(400(7.3263)(4)^2=46,888\), and \(147.1518(7.3263)(3)^2=9,703\), approximately. The best of these candidates is \((0,1)\).

For proportional error, the corresponding determinants are \(102^2(1)^2=10,404\), \(102^2(4)^2=166,464\), and \(102^2(3)^2=93,636\). The best pair is \((0,4)\): with equal weights, the greatest time separation wins.

For combined error, the determinants are approximately \(4.897\times10^3\), \(8.889\times10^3\), and \(3.706\times10^3\), respectively. For example, the \((0,4)\) value is \(81.28(6.8354734)(4)^2=8889.4\), or \(8.889\times10^3\) to four significant figures. The best candidate pair is again \((0,4)\), but its determinant is far below the proportional-error value because the late sample has low concentration and hence low weight under combined error. These are comparisons only among the stated candidate schedules; they do not establish a continuous-time optimum or address operational constraints.

Checks, Failure Modes, and Alternatives

A practical design calculation should be checked at both the model and numerical levels. First verify that each variance is positive and has the appropriate units: \(v\) is in concentration squared, while \(\sigma_p\) is dimensionless. Then verify the sensitivity derivatives and whether the variance depends on the PK parameters. Finally, calculate determinants using unrounded intermediate values where possible; a rounded table can otherwise obscure close schedule rankings.

  • Dropping the variance term automatically. For additive error with fixed \(\sigma_a\), the variance derivative is zero. For proportional and combined error it is generally not. The full normal-likelihood information includes both terms.
  • Assuming low concentration always means low information. Under proportional error, the relative precision is constant, and the log-amplitude/log-rate weights in this example are constant across time. Under additive or combined error, low concentrations usually have less mean information.
  • Comparing determinants across parameterizations without care. A determinant changes under a parameter transformation. Compare candidate times using the same parameterization, or use a justified invariant criterion as in the earlier material on parameterization invariance.
  • Trusting a design based on a poor residual model. Residual plots may show increasing spread, curvature, outliers, or concentration-dependent bias inconsistent with the assumed normal variance. A schedule ranking can be fragile if these features are ignored.

An alternative is to use only the mean-sensitivity contribution \(\sum_i f_\theta f_\theta^\mathsf{T}/v_i\). This is often used as a weighted sensitivity approximation, and it is exact for the mean component; it omits information available through parameter-dependent variance. A different alternative is to compare designs by the precision of one particular parameter rather than the joint determinant. That answers a different design question. For either approach, the residual variance must still be represented consistently.

The calculation also has limitations. It is local to nominal PK parameters, relies on the chosen residual model, and treats error coefficients as known. If those coefficients are uncertain or estimated jointly, their parameters and cross-information must be included. Very small predicted concentrations can also raise below-quantification and assay issues not represented by a normal continuous observation. A high determinant does not guarantee good precision for every individual parameter or protect against model misspecification.

Worked Example: Quantifying the Cost of Omitting Variance Information

Return to the combined-error example at times 0 and 4 hours. The mean-only weights are \(80.0\) and \(100e^{-4}/(0.25+e^{-4})=6.826154\), so the mean-only determinant is \(80.0(6.826154)(4)^2=8,737.48\). With the full likelihood, the weights are \(81.28\) and \(6.835473\), giving \(81.28(6.835473)(4)^2=8,889.4\).

The variance contribution raises this candidate schedule's determinant by about \(151.9\), or \(151.9/8737.48=1.74\%\). The change is modest here, but it is not identically zero, and its relative size can differ across parameter values and candidate times. A sound comparison states whether it uses full Gaussian information or the mean-only approximation rather than silently switching between them.

Key takeaway: Under a Gaussian residual model, calculate the schedule information from both the mean and variance derivatives whenever variance depends on PK parameters. Then compare determinants for candidate schedules using the same nominal model, error coefficients, sample count, and parameterization.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

  1. Under a fixed additive-error variance, which term in the one-observation information formula is zero, and why?
  2. For \(v=\sigma_a^2+\sigma_p^2f^2\), define \(r=\sigma_p^2f^2/v\). Show that the variance-information contribution to the log-amplitude diagonal element is \(2r^2\).
  3. In the two-parameter exponential example, why can the proportional-error model favor widely separated times even when the later concentration is much smaller?
  4. For two samples with weights 5 and 2 at times 1 and 3, calculate the determinant of the information matrix.
  5. Name two assumptions or limitations that should be checked before interpreting a local D-optimal schedule.

Solutions. (1) The term \(v_jv_k/(2v^2)\) is zero because \(v=\sigma_a^2\) is fixed, so every \(v_j=0\). (2) Since \(v_{\log A}=2\sigma_p^2f^2\), the variance term for this diagonal element is \(\tfrac12(2\sigma_p^2f^2/v)^2=2r^2\). (3) With proportional error, \(f^2/v=1/\sigma_p^2\) and \(r=1\), so the weights are constant in time for these parameters; the determinant is therefore proportional to squared time separation. (4) The determinant is \(5(2)(3-1)^2=40\). (5) Examples include the assumed Gaussian, independent residuals; positive correctly specified variance; fixed known error coefficients; and dependence on local nominal PK parameters. A schedule can be optimal under these assumptions yet perform poorly if they are wrong.

Clinical Biostats