Worked HMC Analysis of an Oral PK Model

From an Oral PK Model to a Bayesian HMC Fit

A Bayesian HMC analysis links three pieces: the concentration-time model, the population distribution of individual parameters, and the observed-data likelihood. For an oral dose, clearance, volume, and absorption all influence the concentration curve. Their individual values may also vary together across people. HMC can sample the resulting joint posterior, but a set of numerical estimates is informative only if the model is stated clearly and the chains explore that posterior adequately.

The worked analysis below is entirely hypothetical. It illustrates a single-dose, one-compartment model with first-order absorption and elimination, known dose, and complete bioavailability fixed at \(F=1\). The example is deliberately small enough to follow by hand; it is not a substitute for fitting and checking a real data set.

The expression follows from the standard amount balance for first-order absorption: the gut amount declines at rate \(k_{a,i}\), while the central amount receives absorbed drug and loses drug by elimination at rate \(k_i\). With \(A_g(0)=D_i\) and central amount \(A_c(0)=0\), solving these linear differential equations gives \(A_c(t)=D_i k_{a,i}\{\exp(-k_i t)-\exp(-k_{a,i}t)\}/(k_{a,i}-k_i)\). Dividing by \(V_i\) gives the stated concentration. The limiting expression for equal rates follows by taking the limit as \(k_{a,i}\) approaches \(k_i\), or by applying the derivative of the exponential term.

Population Effects, Correlation, and Observation Model

A convenient population model places normal distributions on the logarithms of the three positive individual parameters. This guarantees positive individual clearance, volume, and absorption rates. Let \(\eta_i=(\eta_{CL,i},\eta_{V,i},\eta_{ka,i})^\mathsf{T}\) be a subject’s random-effect vector. Its components are allowed to be correlated.

Here \(\mu_{CL},\mu_V,\mu_{ka}\) are population means on the log scale; exponentiating one gives the population median (often called the typical value) of that parameter under this model. The vector \(\omega\) contains the three random-effect standard deviations, and \(R\) is their correlation matrix. Thus, for example, \(R_{CL,V}<0\) means subjects with above-typical log clearance tend, in this model, to have below-typical log volume. It does not establish a causal relationship.

For positive concentration measurements, suppose the residual error is proportional on the original scale, or equivalently normal and additive on the log scale:

The residuals are assumed conditionally independent across observations given the individual effects and model parameters. This model presumes positive measured concentrations and does not include a separate additive-error component or below-quantification observations. A Bayesian fit combines this likelihood with proper priors for the population parameters, random-effect scales, correlation matrix, and residual scale. In the example, the priors are centered on plausible positive PK values; random-effect standard deviations and \(\sigma\) have half-normal priors, and \(R\) has an LKJ prior favoring but not forcing moderate correlations. Priors must be specified on the scale actually used and checked for plausible implications for predicted concentrations.

Likelihood Contribution and HMC Sampling

For a given set of parameters and random effects, the likelihood compares each observed concentration with its model prediction. Under the lognormal observation model, the log-scale residual is \(r_{ij}=\log(Y_{ij}/C_{ij})\). Its conditional distribution is normal with mean zero and standard deviation \(\sigma\). Because the data are modeled as concentrations rather than log concentrations, the density on the concentration scale also includes the lognormal Jacobian factor \(1/Y_{ij}\); omitting it would change the likelihood used for inference.

HMC samples the joint posterior over population parameters, random effects, and error parameters by using gradients of the log posterior. In practice, positive parameters are represented through unconstrained transformations, and the correlation matrix is parameterized so that sampled values remain valid correlation matrices. The sampling target is the posterior implied by the specified likelihood and priors; HMC does not verify that the PK model or priors are scientifically adequate.

For a worked hypothetical data set, suppose six subjects each receive a single 400 mg oral dose and have concentrations measured at 0.5, 1, 2, 4, and 8 hours. Assume the times and doses are known, all recorded concentrations are positive, and there are no missing or below-quantification values. A fit might use four chains, each with 1,000 warmup iterations followed by 2,000 retained iterations. Warmup adapts the numerical integration settings; those draws are not included in posterior summaries.

Summarizing Population and Correlation Results

Assume the hypothetical run has posterior medians and central 90% credible intervals as shown below. Each interval is a posterior interval conditional on the model, priors, and data. It is not a statement that the parameter has a 90% frequentist coverage probability, nor is it a probability that a fixed parameter changes value.

The correlation interval spans zero, so this hypothetical fit does not sharply determine the direction of the clearance-volume random-effect association. A negative posterior median alone is not strong evidence that the population correlation is negative. The uncertainty in a correlation can be substantial, especially when the number of subjects is small or the individual parameters are weakly identified by the sampling times.

Diagnostics, Alternatives, and Limitations

A useful review separates sampling diagnostics from model checks. HMC diagnostics ask whether the algorithm appears to explore the target posterior. Model checks ask whether that posterior model can reproduce relevant features of the data. For the hypothetical oral PK fit, inspect trace plots and convergence summaries, then compare observed concentrations with posterior predictive concentrations across subjects and time. Residual patterns that systematically rise or fall with time may indicate a misspecified structural or residual model even when chains mix well.

The model also has substantive limitations. Fixing \(F=1\) treats bioavailability as known; with oral data alone, \(F\) and \(V\) can be difficult to distinguish because concentration depends on their ratio. A single dose and sparse sampling may not separately identify absorption, elimination, and variability. The one-compartment, first-order absorption structure excludes distribution phases, lag time, time-varying clearance, and more complex residual error. The posterior is conditional on these simplifications and on the priors.

Bayesian HMC provides joint posterior distributions and directly represents uncertainty and dependence among parameters. Point-estimation alternatives such as likelihood-based FOCE or SAEM can be computationally efficient and familiar in population PK workflows, but typically report estimates and approximations to uncertainty rather than the full Bayesian posterior. Neither approach rescues a weakly informative design or a poor structural model; estimates from either should be checked against the data and scientific assumptions.