MCMC for Individual PK Posterior Sampling

Why Sample an Individual PK Posterior?

An individual concentration record may constrain PK parameters only imperfectly. A nonlinear concentration model can also make the posterior asymmetric or otherwise poorly represented by a single point estimate and a local approximation. Markov chain Monte Carlo (MCMC) addresses this by generating parameter values whose long-run distribution is the individual posterior. The resulting draws can summarize uncertainty in clearance or support later forecasting; they do not make weakly informative data more informative.

This tutorial develops a Metropolis–Hastings sampler for a simple individual PK model, then shows how to assess its output. The same logic extends to a vector of individual parameters, although dependence among parameters can make sampling substantially more difficult.

Target Posterior and Metropolis–Hastings Sampling

Let \(k=CL/V\) be an elimination rate constant in \(\mathrm{h}^{-1}\). For an illustrative one-compartment IV bolus model with known volume \(V\) and dose \(D\), the predicted concentration at time \(t\) is \(C(t)= (D/V)\exp(-kt)\). For observations \(y_i\) at times \(t_i\), suppose, for this example, that residuals are independent Normal with known standard deviation \(\sigma\). With prior density \(p(k)\), the posterior density is proportional to the likelihood times the prior:

The proportionality omits constants that do not depend on \(k\). MCMC needs the posterior density only up to such a constant. Since \(k\) must be positive, define \(x=\log k\), so every real value of \(x\) maps to a valid \(k\). The target for a chain on \(x\) must be the posterior density with respect to \(x\). If starting from a density with respect to \(k\), this target includes the Jacobian factor \(dk/dx=\exp(x)\). Equivalently, specify a prior directly on \(x\), as in the worked example below.

A common choice is \(x'=x+\epsilon\), where \(\epsilon\sim N(0,s_q^2)\) and \(s_q\) is a proposal scale. The transition rule allows moves to lower-density regions, which is essential: repeatedly moving only uphill would not generate the full posterior. A proposal’s acceptance probability is not itself the probability that the current state is correct; it is part of a rule that produces the desired stationary distribution.

Why the Transition Rule Targets the Posterior

For a symmetric proposal, the acceptance probability is \(\min\{1,\pi(x')/\pi(x)\}\). The probability flow from \(x\) to \(x'\), when weighted by the target density, is proportional to \(\pi(x)q(x'\mid x)\min\{1,\pi(x')/\pi(x)\}=q(x'\mid x)\min\{\pi(x),\pi(x')\}\). By proposal symmetry this equals the flow from \(x'\) to \(x\). Thus the transition satisfies detailed balance with respect to \(\pi\), which makes \(\pi\) stationary. This argument establishes stationarity, not that a finite run has reached it. A chain must also be able to explore the relevant target support; poor proposals, constraints, or separated regions can prevent useful exploration in practice.

For a vector of log-transformed PK parameters, the same acceptance rule applies to the joint posterior density. A proposal may change one parameter at a time or propose a vector move. The latter can be more efficient when parameters are strongly dependent, but selecting an effective proposal can be challenging. This tutorial’s scalar example isolates the sampling logic rather than resolving the broader issue of posterior dependence.

Burn-In, Mixing, and Effective Sample Size

A chain started at an arbitrary point may initially reflect its starting value more than the target posterior. Burn-in is an initial segment that is excluded from summaries to reduce this initialization effect. Choosing a burn-in length does not prove convergence: a chain can appear stable while remaining trapped in a region it has not adequately explored. Multiple chains started at dispersed plausible values help reveal failures that a single trace may conceal.

Mixing describes how rapidly a chain moves through the target distribution. Strong autocorrelation means nearby retained states carry overlapping information. If \(M\) retained draws have lag-\(j\) autocorrelation \(\rho_j\), a common approximation to effective sample size is \(ESS\approx M/[1+2\sum_{j=1}^{\infty}\rho_j]\), using a suitable truncation of the estimated autocorrelation sequence in practice. ESS is an approximate count of independent draws with comparable precision for a particular summary; it is not the number of distinct parameter values and can differ for different summaries.

Diagnostics, Alternatives, and Limitations

A trace plot should show the sampled parameter against iteration. Warning signs include persistent drift, long nearly constant stretches, chains occupying different ranges, or abrupt behavior suggesting numerical or implementation errors. Autocorrelation plots help identify slow mixing; a low acceptance rate can signal proposals that are too large, while an extremely high rate can signal tiny steps and slow exploration. Neither acceptance rate alone is a convergence test. Repeatedly evaluating the likelihood at known parameter values and checking the transformed-scale Jacobian are useful implementation checks.

MCMC is useful when posterior summaries are needed and a nonlinear posterior is not adequately captured by a local approximation. A MAP estimate gives a posterior mode, and a Laplace approximation summarizes the posterior locally around that mode. Those methods can be quicker and useful when the approximation is suitable; MCMC instead represents the target through a sample, at the cost of computation and the need to assess convergence and Monte Carlo error. MCMC does not automatically solve weak identification, prior-data conflict, or model misspecification. Its sample can faithfully represent a posterior that is broad or sensitive to assumptions.

A practical workflow is to define the posterior and its parameter scale, validate the target calculation, run multiple chains, inspect traces and diagnostics, and increase simulation effort when mixing or Monte Carlo precision is inadequate. Only after those checks should retained draws be used to report individual parameter summaries or downstream predictions. For a dosing decision, report the clinically relevant posterior quantity and its units, not just an iteration count or an acceptance rate.

Common Pitfall

• Wrong reasoning: “The chain accepted many proposals and its trace looks roughly stable, so the posterior has converged.”
• Correction: Acceptance and a visually stable trace are limited evidence. Compare multiple dispersed chains, check autocorrelation and relevant convergence diagnostics, and assess effective sample size and Monte Carlo error for the actual summaries being used.
• Wrong reasoning: “Burn-in removes the effect of initialization, so all draws after it are independent.”
• Correction: Burn-in may reduce initialization effects, but retained draws can remain highly autocorrelated. Their effective sample size can be far smaller than their count.