Tutorials › Optimal Design and Model-Informed Dosing › Numerical Search Over Candidate Sampling Times

Optimal Sampling for Individual PK Models · Tutorial 48 of 497

Numerical Search Over Candidate Sampling Times

Implement and verify numerical searches for D-optimal sampling times, while checking whether a solution is stable, nonsingular, and plausibly optimal.

Advanced 12 min read

What You'll Learn

  • Construct a local Fisher information matrix from model sensitivities at proposed sampling times.
  • Compare fixed-size schedules by their log-determinants without enumerating large determinants directly.
  • Use exhaustive candidate searches or exchange algorithms, depending on the size of the candidate set.
  • Search over continuous sampling times from multiple starting schedules and compare the resulting objectives.
  • Diagnose singular information matrices, premature stopping, and apparent local optima.
  • Interpret the limitations of a numerical optimum as evidence about the specified model and design problem.

Why Numerical Search Needs Diagnostics

For an individual PK design, a sampling schedule can be represented by its collection of observation times. The Fisher information at those times depends on how strongly and distinctly the predicted observations respond to the parameters. D-optimality ranks schedules by the determinant of that information matrix: larger determinants indicate greater joint information about the parameter vector, under the model and assumptions used to calculate the matrix.

In a small candidate set, it may be possible to calculate the determinant for every schedule. With many candidate times, or when times can vary continuously, numerical optimization is more practical. But an optimizer returning a schedule is not proof that it found a global optimum. A sound search therefore combines a stable objective calculation with multiple starting schedules, checks for numerical rank, and stopping diagnostics that examine both the objective and the search itself.

Definition: For a proposed schedule \(\xi\), let \(\mathcal I(\xi)\) be its Fisher information matrix for the parameter vector of interest. A nonsingular schedule is D-optimal over a specified design space if it maximizes \(\det\{\mathcal I(\xi)\}\) over that space. Numerical searches commonly maximize \(\log\det\{\mathcal I(\xi)\}\), which has the same maximizer when the information matrix is positive definite.

The search is conditional on the model, nominal parameter values when the design is local, the observation-error assumptions, and the allowed set of schedules. It does not establish that those inputs are correct. In the example below, the model is deliberately simple so the numerical calculations can be checked by hand.

Information and a Stable Objective

Suppose log concentration at time \(t\), measured in hours, follows the one-compartment IV bolus form \(Y(t)=a-kt+\epsilon\), where \(a\) is log concentration at time zero, \(k\) is the elimination rate constant in \(\mathrm{h}^{-1}\), and the independent errors satisfy \(\epsilon\sim N(0,\sigma^2)\). Assume that \(\sigma=0.1\) log-concentration units is known and common to all observations. For a fixed schedule, the sensitivity of the mean response to \((a,k)\) is \((1,-t)\). Under these assumptions, the information matrix is the sum of the sensitivity outer products, scaled by \(1/\sigma^2\).

$$ \mathcal I(\xi) = \frac{1}{\sigma^2} \sum_{i=1}^{n} \begin{pmatrix}1\\-t_i\end{pmatrix} \begin{pmatrix}1&-t_i\end{pmatrix} = \frac{1}{\sigma^2} \begin{pmatrix} n&-\sum_i t_i\\ -\sum_i t_i&\sum_i t_i^2 \end{pmatrix}. $$

The matrix is positive definite exactly when the schedule contains at least two distinct times. Its determinant can be written \(\{n\sum_i t_i^2-(\sum_i t_i)^2\}/\sigma^4\). This identity both gives a useful hand-check and exposes a failure mode: repeated identical times make the sensitivity rows identical, so the information matrix is singular.

For a general positive-definite information matrix, a numerical implementation should usually obtain the log-determinant through a Cholesky factorization, \(\mathcal I=LL^\mathsf{T}\). Then \(\log\det(\mathcal I)=2\sum_j\log L_{jj}\). This avoids forming a potentially enormous or extremely small determinant. If Cholesky factorization fails, the schedule may be singular or nearly singular; silently adding a small number to the diagonal can hide that design problem and should not be treated as a valid fix.

Searching a Finite Set of Candidate Times

When the candidate set is modest, exhaustive search is transparent: list every admissible fixed-size schedule, calculate its information matrix, and compare log-determinants. If there are \(m\) candidate times and \(n\) distinct sampling times are required, there are \(\binom{m}{n}\) schedules. For a larger set, an exchange search can begin with one schedule, replace one selected time by an unselected candidate when that improves the objective, and repeat until no tested single replacement improves it. Such a search is efficient but may stop at a schedule that is only locally best under the allowed replacements.

Worked Example: Exhaustive Search on a Four-Time Grid

Consider three distinct samples chosen from the hypothetical candidate times \(0,1,2,3\) hours under the model above. The assumptions are independent normal errors with known \(\sigma=0.1\), and the two parameters are \(a\) and \(k\). For \(n=3\), the determinant is \( \{3\sum t_i^2-(\sum t_i)^2\}/0.1^4\).

For schedule \((0,1,2)\), \(\sum t_i=3\) and \(\sum t_i^2=5\). Thus the numerator is \(3(5)-3^2=6\), and \(\det(\mathcal I)=6/0.0001=60{,}000\). Its log-determinant is \(\log(60{,}000)\approx11.0021\).

For schedule \((0,1,3)\), the corresponding sums are \(4\) and \(10\). The numerator is \(3(10)-4^2=14\), so the determinant is \(14/0.0001=140{,}000\) and the log-determinant is approximately \(11.8494\). The schedules \((0,2,3)\) and \((1,2,3)\) give numerators \(14\) and \(6\), respectively. Therefore the two tied best schedules are \((0,1,3)\) and \((0,2,3)\).

The conclusion follows from checking all four possible schedules, not merely from observing that one selected schedule improved. Relative to a determinant of \(140{,}000\), a determinant of \(60{,}000\) has D-efficiency \((60{,}000/140{,}000)^{1/2}=\sqrt{3/7}\approx0.6547\) for this two-parameter problem. The candidate-grid result is optimal only among these four schedules; it does not by itself establish the best schedule among all real-valued times.

Searching Continuous Sampling Times

For continuous search, the optimizer varies time values within the permitted time domain rather than selecting only from a grid. The parameterization should ensure that the proposed times remain in the domain and, when relevant, in the intended order. One can use a bounded optimizer or transform unconstrained variables into valid times. The transformed objective must still be evaluated from the actual times and information matrix, not from rounded display values.

The following calculation illustrates the objective and its boundary behavior. With two observations at \(t_1<t_2\), the determinant above simplifies to \((t_2-t_1)^2/\sigma^4\). Consequently the log-determinant is \(2\log(t_2-t_1)-4\log(\sigma)\). This is increasing in the time separation. Its derivatives with respect to \(t_1\) and \(t_2\) are \(-2/(t_2-t_1)\) and \(2/(t_2-t_1)\), respectively. Thus, on a bounded domain, the objective drives the earlier time toward the lower endpoint and the later time toward the upper endpoint; an interior zero-gradient solution is not expected.

Worked Example: Multiple Starts on a Continuous Interval

Suppose two samples are to be placed in the interval from 0 to 4 hours, with \(t_1<t_2\). The same hypothetical log-concentration model, known error standard deviation \(\sigma=0.1\), and independent-error assumptions apply. Start a bounded search from \((0.5,2.0)\), \((1.0,3.5)\), and \((0.2,3.8)\) hours.

For the first start, the gap is \(1.5\) hours, so \(\det(\mathcal I)=1.5^2/0.1^4=22{,}500\) and \(\log\det(\mathcal I)\approx10.0213\). The second start has gap \(2.5\), giving determinant \(62{,}500\) and log-determinant approximately \(11.0429\). The third has gap \(3.6\), giving determinant \(129{,}600\) and log-determinant approximately \(11.7722\).

A bound-aware search moves the earlier time toward 0 and the later time toward 4. At \((0,4)\), the gap is 4, so the determinant is \(4^2/0.1^4=160{,}000\), and the log-determinant is \(\log(160{,}000)\approx11.9829\). Since the derived objective increases with the gap and no gap can exceed 4 on this interval, this result is globally optimal for the stated problem. In a more complicated PK model, agreement across starts is useful evidence of stability but does not provide this kind of proof.

Stopping, Rank, and Failure Diagnostics

A numerical stopping rule should be based on objective changes at adequate precision, not on a rounded determinant in a report. For a maximization algorithm, a practical check is that the log-determinant improvement over successive iterations is below a chosen tolerance and that the design variables or proposed exchanges have also stopped changing materially. In continuous optimization, inspect a suitable gradient or projected-gradient condition as well: at an interior optimum, the gradient should be near zero; at a boundary optimum, the gradient need not be zero, and its direction must be consistent with the active bounds.

Worked Example: Detecting a Singular Schedule

Take two observations at the same time, \(t_1=t_2=1\) hour. The sensitivity rows are both \((1,-1)\), so the information matrix is \(\mathcal I=100\begin{pmatrix}2&-2\\-2&2\end{pmatrix}\). Its determinant is \(100^2(4-4)=0\). There is no finite log-determinant, and the two parameters cannot be separately informed by these identical sensitivity directions.

As a check, change the schedule to \((1,2)\) hours. The determinant formula gives \((2-1)^2/0.1^4=10{,}000\), with log-determinant \(\log(10{,}000)\approx9.2103\). The second schedule is nonsingular because the two sensitivity rows differ. This calculation does not say that the information is adequate for every practical purpose; it only confirms that the matrix has full rank in this example.

Useful diagnostics include verifying that the information matrix is symmetric, checking its numerical rank or Cholesky factorization, and confirming that all reported times satisfy the design rules supplied to the optimizer. Compare high-precision objective values across starts and, where possible, against an independent calculation or a finer candidate grid. If several starts stop at different log-determinants, investigate whether the optimizer has found different local optima, reached different boundaries, or failed because of poor scaling. If all starts give the same rounded determinant but differ at meaningful precision, the apparent agreement may be an artifact of rounding.

A finite-grid search, an exchange algorithm, and continuous optimization answer related but distinct computational questions. Exhaustive search guarantees the best schedule on a fully enumerated finite set, but its cost grows combinatorially. Exchange methods can be faster on large grids, but their stopping condition usually certifies only that no tested local exchange improves the objective. Continuous methods can explore between grid points, but may depend on starting values and optimizer settings. A fine grid is a useful cross-check, not a proof that a continuous optimum has been found. Approximate designs, which allow support weights rather than requiring a fixed list of individual observations, are a different design representation and should not be confused with an exact fixed-size schedule.

Common Pitfalls and Interpretation

  • Using the determinant after rounding. Calculate the objective from full-precision matrix entries, then round the reported result. Early rounding can change a displayed determinant or even the ranking of close schedules.
  • Treating one optimizer run as a global certificate. Run multiple starts, compare objective values, and inspect the final schedules. A local algorithm can stop without finding the global maximum.
  • Ignoring singularity. A determinant of zero means the matrix is rank deficient; it is not an ordinary low but usable log-determinant. Diagnose the sensitivities and schedule rather than concealing the problem by arbitrary numerical regularization.
  • Confusing numerical convergence with model adequacy. A stable maximum is optimal only for the specified model, error structure, parameter values, and design space. It does not validate those assumptions or guarantee good estimation under model misspecification.

The determinant is sensitive to parameter units, while a smooth reparameterization with a nonsingular Jacobian, applied consistently at the same parameter point, multiplies all candidate determinants by the same positive factor and therefore preserves their ranking. Regardless, report parameterization and units clearly when communicating numerical values. The most useful result is not an isolated determinant: it is a schedule, its objective at sufficient precision, the assumptions that generated its information matrix, and diagnostics that show how convincingly the search stopped.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. For the two-parameter log-concentration model, why does a schedule with two identical sampling times have a singular information matrix?
  2. Three samples are selected from candidate times \(0,1,2,3\). What is the determinant for schedule \((0,2,3)\) when \(\sigma=0.1\)?
  3. Why is maximizing \(\log\det(\mathcal I)\) equivalent to maximizing \(\det(\mathcal I)\) for positive-definite information matrices?
  4. What should be checked if different continuous-search starts return different objective values?
  5. Does a fine candidate grid prove that a reported schedule is the continuous-time global optimum? Explain.

Solutions. (1) The sensitivity rows \((1,-t)\) are identical at identical times, so the information matrix has rank one rather than two. (2) For \((0,2,3)\), \(\sum t_i=5\), \(\sum t_i^2=13\), and \(3(13)-5^2=14\); therefore the determinant is \(14/0.0001=140{,}000\). (3) The logarithm is strictly increasing for positive numbers, so it preserves their ordering. (4) Compare the final schedules and high-precision log-determinants, inspect boundaries and stopping diagnostics, and consider local optima, scaling, or convergence failure. (5) No. A grid checks only its discrete candidate schedules; continuous times between grid points may yield a better value, and a continuous-search result itself needs suitable diagnostics or a global justification.

Clinical Biostats