Sequential Bayesian Forecasting with Therapeutic Drug Monitoring

Why Forecasts Are Updated Sequentially

Therapeutic drug monitoring (TDM) adds information about an individual’s pharmacokinetics as treatment proceeds. A forecast made before a concentration is measured reflects the population prior and any earlier measurements. Once a new result is available, the parameter distribution and the forecast for a later concentration should both be revised. The update is not simply a correction to the previous point forecast: the new result changes how plausible different individual parameter values are, and therefore changes the full predictive distribution.

The focus here is the sequence of updates for one individual under a specified regimen and observation model. The population model supplies the prior distribution for that individual’s parameters; the likelihood connects measured concentrations to those parameters. The introductory material on individual parameter updating and posterior predictive distributions establishes these ingredients. Here, the emphasis is on applying them recursively and checking how each new result changes the forecast.

Definition: Let \(\boldsymbol{\theta}\) denote an individual’s pharmacokinetic parameters, \(D_k=(z_1,\ldots,z_k)\) the first \(k\) monitoring results on the analysis scale, and \(x_{1:k}\) their known dosing, sampling-time, and covariate information. The sequential parameter posterior is \(p(\boldsymbol{\theta}\mid D_k,x_{1:k})\). For a future concentration \(C_*\), the posterior predictive distribution is the distribution of \(C_*\) conditional on the observed results and the known future conditions.

The Recursive Update

Suppose the model specifies the conditional density of the next monitoring result given the individual parameters and the relevant dosing and sampling information. If each result is conditionally independent of earlier results given those parameters and the modeled information, Bayes’ rule gives the update

p(\boldsymbol{\theta}\mid D_k,x_{1:k}) \propto p(z_k\mid\boldsymbol{\theta},x_k)\, p(\boldsymbol{\theta}\mid D_{k-1},x_{1:k-1}).

The previous posterior is the next update’s prior. The proportionality constant is chosen so the updated density integrates to one. In a nonlinear model, \(p(z_k\mid\boldsymbol{\theta},x_k)\) is obtained from the model-predicted concentration and the observation-error distribution; the update will not generally have a closed-form solution. Numerical integration, sampling, or another suitable approximation may be used to represent the posterior.

For a future concentration, parameter uncertainty and future observation variability both matter. If \(f(C_*\mid\boldsymbol{\theta},x_*)\) is the model’s conditional predictive density at the planned future conditions, then

p(C_*\mid D_k,x_{1:k},x_*) = \int f(C_*\mid\boldsymbol{\theta},x_*)\, p(\boldsymbol{\theta}\mid D_k,x_{1:k})\,d\boldsymbol{\theta}.

Thus, each new TDM result first updates the parameter posterior and then changes the distribution being averaged over to forecast \(C_*\). The predictive distribution also includes any future residual variation represented by \(f\). This is why the parameter posterior alone is not a prediction interval for a future concentration.

Assumptions for the recursive update: The population prior and concentration model are appropriate for the individual; dosing, sampling times, and relevant covariates are recorded correctly; the observation-error model adequately represents measurement and residual variability; and the conditional likelihood accounts for dependence among results. The simple product of one-result likelihoods is valid only under conditional independence given the modeled parameters and information. If residuals are correlated, the joint likelihood must represent that dependence instead.

A Gaussian Benchmark for the Update

A conjugate calculation makes the sequential logic transparent. Let \(\theta\) be a scalar individual parameter and suppose the transformed result \(y_i\) has the model

y_i=\theta+\epsilon_i,\qquad \theta\sim N(m_0,v_0),\qquad \epsilon_i\sim N(0,s_e^2),

where the errors are independent of one another and of \(\theta\), with known variance \(s_e^2\). After \(k\) results, the posterior is normal with precision equal to prior precision plus the information from the \(k\) observations:

v_k^{-1}=v_0^{-1}+\frac{k}{s_e^2}, \qquad m_k=v_k\left(\frac{m_0}{v_0}+\frac{\sum_{i=1}^k y_i}{s_e^2}\right).

To justify the result, multiply the prior density by the observation likelihoods. Ignoring terms that do not depend on \(\theta\), the exponent is

-\frac{1}{2}\left[ \frac{(\theta-m_0)^2}{v_0} +\sum_{i=1}^{k}\frac{(y_i-\theta)^2}{s_e^2} \right].

Collecting the terms involving \(\theta^2\) gives coefficient \(v_0^{-1}+k/s_e^2\); collecting the terms involving \(\theta\) gives \(m_0/v_0+\sum_i y_i/s_e^2\). Completing the square therefore produces the stated variance and mean. This benchmark requires a linear observation relation on the analysis scale and known variances. It is a useful check on sequential calculations, not a general formula for nonlinear pharmacokinetic models.

Worked Example: Updating After the First Result

Consider a hypothetical individual whose log-clearance parameter has prior distribution \(\theta\sim N(0.50,0.09)\). For illustration, use a simplified log-concentration model \(z_i=2-\theta+\epsilon_i\), where \(\epsilon_i\sim N(0,0.04)\). The known dose and sampling-time effects are represented by the offset 2 in this numerical example. The errors are assumed independent, and all log concentrations are defined as \(z_i=\ln(C_i/(1\ \mathrm{mg/L}))\), so their logarithms are dimensionless.

For the first result \(z_1=1.60\), rearrange the observation model as \(y_1=2-z_1=\theta+\epsilon'_1\), where \(\epsilon'_1=-\epsilon_1\) still has variance \(0.04\). Thus \(y_1=0.40\). The prior precision is \(1/0.09=11.111\), and the observation precision is \(1/0.04=25\). The posterior variance and mean are

v_1=\frac{1}{11.111+25}=0.02769,\qquad m_1=\frac{11.111(0.50)+25(0.40)}{36.111}=0.4308.

After this result, the individual’s posterior is approximately \(N(0.4308,0.02769)\). The observed concentration is \(\exp(1.60)=4.953\ \mathrm{mg/L}\), but that measurement is not itself the forecast for a later sample: the forecast requires propagating the updated parameter distribution through the future concentration model and including future residual variability.

Following Both the Point Forecast and Its Interval

Continue with the same hypothetical model and suppose a second result is \(z_2=1.30\), giving \(y_2=2-1.30=0.70\). The update can be made from the first posterior: add 25 to its precision and add \(25(0.70)\) to its precision-weighted mean. Equivalently, use both observations in the batch formula. The result is

v_2=\frac{1}{11.111+25+25}=0.01636,\qquad m_2=\frac{11.111(0.50)+25(0.40)+25(0.70)}{61.111}=0.5409.

At a future condition represented by the same offset of 2, the future log concentration is \(z_*=2-\theta+\epsilon_*\), with independent \(\epsilon_*\sim N(0,0.04)\). Integrating over the normal posterior gives a normal predictive distribution with mean \(2-m_k\) and variance \(v_k+0.04\). The \(v_k\) term is parameter uncertainty; the \(0.04\) term is future residual variability.

Worked Example: Revising the Forecast After Two Results

Using the posterior after the first result, the future log-concentration forecast has mean \(2-0.4308=1.5692\) and variance \(0.02769+0.04=0.06769\), giving a standard deviation of \(0.2602\). Its 95% predictive interval on the log scale is \(1.5692\pm1.96(0.2602)\), or approximately \([1.059,2.079]\).

After the second result, the predictive mean becomes \(2-0.5409=1.4591\), and its variance is \(0.01636+0.04=0.05636\), giving a standard deviation of \(0.2374\). The 95% predictive interval on the log scale is

1.4591\pm1.96(0.2374)=[0.9938,1.9244].

Exponentiating the endpoints gives a 95% predictive interval of approximately \([2.701,6.851]\ \mathrm{mg/L}\). This interval is for a future measured concentration under the stated model and conditions. It is not an interval for the parameter \(\theta\).

For a point forecast on the concentration scale, use the mean of the lognormal predictive distribution, not simply the exponential of its log-scale mean. If \(Z\sim N(\mu,v)\), then \(E[\exp(Z)]=\exp(\mu+v/2)\). After the second result, the forecast mean concentration is therefore \(\exp(1.4591+0.05636/2)=4.425\ \mathrm{mg/L}\). For comparison, after the first result it was \(\exp(1.5692+0.06769/2)=4.968\ \mathrm{mg/L}\). In this sequence, the second result shifts the point forecast downward and narrows the interval. Neither a shift nor a narrowing is guaranteed for every new observation.

Incorporating a Further TDM Result

The next example continues the same hypothetical sequence and emphasizes the effect of a new result that is relatively high compared with the second measurement. For clarity, retain the simplified model and its assumptions; the arithmetic does not imply that this model is adequate for every drug or sampling schedule.

Worked Example: Updating with a New Concentration of 5.280 mg/L

Suppose the next log result is \(z_3=1.6639\), corresponding to a measured concentration of approximately \(\exp(1.6639)=5.280\ \mathrm{mg/L}\). On the model’s transformed parameter scale, this observation gives \(y_3=2-z_3=0.3361\). Start from the two-result posterior, with mean \(0.5409\) and variance \(0.01636\). Its precision is approximately \(61.111\); adding the new observation precision of 25 gives a posterior precision of \(86.111\).

The updated variance and mean are

v_3=\frac{1}{86.111}=0.01161,\qquad m_3=\frac{61.111(0.5409)+25(0.3361)}{86.111}=0.4814.

The future log-concentration forecast now has mean \(2-0.4814=1.5186\), variance \(0.01161+0.04=0.05161\), and standard deviation \(0.2272\). Its 95% predictive interval on the log scale is \(1.5186\pm1.96(0.2272)=[1.073,1.964]\), approximately. On the concentration scale this is about \([2.925,7.127]\ \mathrm{mg/L}\). The predictive mean concentration is \(\exp(1.5186+0.05161/2)=4.685\ \mathrm{mg/L}\).

Compared with the two-result forecast, this new, higher concentration shifts the forecast mean upward, from \(4.425\) to \(4.685\ \mathrm{mg/L}\). The interval also changes: its endpoints are now approximately \(2.925\) and \(7.127\ \mathrm{mg/L}\). The direction and size of an update depend on the observation, the model, and the previous posterior; a new result does not mechanically pull every forecast toward itself by the same amount.

Diagnostics, Alternatives, and Limitations

A sequential calculation is only as credible as its likelihood and conditioning information. Review the observed result against the model-predicted concentration at that result’s sampling time and dose history. Check units, dose records, sample timing, assay limits, and whether the error scale is appropriate. Large, repeated discrepancies or residual patterns across sampling times can signal a poor structural or error model. A single surprising observation should prompt investigation, not automatic deletion or an unexamined parameter shift.

Also check the posterior itself. In nonlinear models it may be skewed, multimodal, or concentrated near parameter boundaries. A single mean and variance may hide these features. Numerical approximations should be checked for adequate coverage of relevant parameter values and stable predictive summaries. If measurements share residual dependence, treating them as independent can make the posterior too concentrated; the likelihood must represent the dependence or other shared sources of variation.

Sequential updating and batch updating are equivalent when they use the same prior and correctly specified joint likelihood: multiplying the successive likelihood contributions gives the same posterior as using them together. Sequential calculation is convenient when results arrive over time, while batch calculation can be simpler for a fixed data set. A plug-in forecast based on one parameter estimate is easier to compute, but it omits parameter uncertainty and can understate predictive uncertainty, especially when individual data are sparse. Neither an exact recursive formula nor a sophisticated numerical method protects against incorrect data, a misspecified model, or unmodeled changes in the individual.

A further limitation is that intervals are conditional on the model and the information supplied. In the example, the predictive variance decreases after each result because the posterior parameter variance decreases, while the residual variance remains fixed. In general, new observations may increase or decrease predictive uncertainty; and uncertainty not represented in the model will not be captured by the resulting interval. Interpret a nominal 95% predictive interval as a model-based interval intended to contain the future concentration with probability 0.95 under the specified conditioning information and assumptions.

Common pitfall: Treating the latest measured concentration as the next point forecast, or treating the updated parameter posterior as the predictive distribution. The correction is to update the parameter posterior with the new result, then integrate the future concentration model over that posterior and include future residual variability. The resulting predictive forecast may differ from the last observed value and has its own uncertainty.

Check Your Understanding

Answer each question using the definitions and model distinctions developed above. Solutions follow the questions.

Under what condition can the likelihood contributions from successive TDM results be multiplied as separate terms?
In the Gaussian benchmark, what happens to the posterior precision after one more independent result with error variance \(s_e^2\)?
Why does the future predictive variance in the numerical example include both \(v_k\) and \(0.04\)?
After the third result in the worked sequence, is the 95% interval for \(\theta\) or for a future concentration? Give its approximate concentration-scale endpoints.
Why can a new TDM result make a predictive interval wider rather than narrower?

Solutions. (1) Separate likelihood terms are justified when results are conditionally independent given the modeled individual parameters and relevant dosing, sampling, and covariate information. (2) The precision increases by \(1/s_e^2\), because the new independent observation contributes that amount of precision. (3) \(v_k\) represents remaining uncertainty about the individual parameter, while \(0.04\) represents residual variability in the future measurement. (4) It is a future-concentration predictive interval, approximately \([2.925,7.127]\ \mathrm{mg/L}\), not an interval for \(\theta\). (5) A new observation can shift posterior weight toward parameter values whose future predictions are more variable or farther apart; the total predictive uncertainty need not decrease for every observation.