Updating Individual Parameter Distributions with Observed Concentrations

Why an Individual Distribution Changes When Data Arrive

Before an individual has concentration measurements, the population model supplies a distribution of plausible pharmacokinetic parameters for that individual, conditional on the population information and any relevant covariates. Once concentrations are observed, some parameter values predict those data better than others. Updating means combining this population distribution with the likelihood of the individual’s observed concentrations. The result is a distribution, not automatically a single “best” parameter value.

The update is conditional on the model, dose history, sampling times, covariates, and population parameters being used. It does not by itself account for uncertainty in estimated population parameters if those parameters have simply been fixed at fitted values.

Definition: Let $\theta_i$ be the vector of pharmacokinetic parameters for individual $i$, $\psi$ the population-model parameters, and $y_i=(y_{i1},\ldots,y_{im})$ the individual’s observed concentrations at known sampling times and doses. The individual posterior distribution is $$ p(\theta_i\mid y_i,\psi)= \frac{p(y_i\mid\theta_i,\psi)\,p(\theta_i\mid\psi)} {\int p(y_i\mid\vartheta,\psi)\,p(\vartheta\mid\psi)\,d\vartheta}. $$ The prior $p(\theta_i\mid\psi)$ is the population distribution for an individual with the specified covariates; $p(y_i\mid\theta_i,\psi)$ is the concentration likelihood. The denominator normalizes the product to integrate to one.

Likelihood and Assumptions

Write the structural concentration predicted at sample $j$ as $f(t_{ij};\theta_i)$, with the dose history and covariates understood. Under independent additive Gaussian residual error with known standard deviation $\sigma$, the observation model is

y_{ij}=f(t_{ij};\theta_i)+\epsilon_{ij},\qquad \epsilon_{ij}\sim N(0,\sigma^2).

Conditional on $\theta_i$, the observations are assumed independent in this formulation. The likelihood is therefore

p(y_i\mid\theta_i,\psi)= \prod_{j=1}^{m}\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{[y_{ij}-f(t_{ij};\theta_i)]^2}{2\sigma^2}\right\}.

These assumptions must match the data and model. For proportional error, a Gaussian model on log concentrations may be more suitable than additive Gaussian error on the original scale. If residuals are correlated, censored, or have time-varying variance, the likelihood must represent that structure; multiplying independent constant-variance Gaussian terms would then be incorrect. The dose and sampling records must also be correctly specified, since a wrong input history can look like an unusual individual parameter.

When $f$ is nonlinear in $\theta_i$, the posterior is generally not a familiar named distribution. It may be skewed, multimodal, or constrained by parameter bounds. The Bayes formula still applies, but its integral and summaries may require numerical methods. The next examples first show exact Gaussian updates, then a nonlinear likelihood calculation.

A Gaussian Update and Its Derivation

A useful closed-form case arises when a scalar individual parameter $\eta_i$ has a normal population prior and the data reduce to Gaussian observations that are linear in $\eta_i$. Suppose $\eta_i\sim N(\mu_0,\tau^2)$ and $z_j=\eta_i+\epsilon_j$, where the independent errors have known variances $\sigma_j^2$. The prior precision is $1/\tau^2$; each observation contributes precision $1/\sigma_j^2$.

p(\eta_i\mid z_1,\ldots,z_m)\ \propto\ \exp\left\{-\frac12\left[ \frac{(\eta_i-\mu_0)^2}{\tau^2} +\sum_{j=1}^m\frac{(z_j-\eta_i)^2}{\sigma_j^2} \right]\right\}.

To identify the posterior, collect the terms involving $\eta_i^2$ and $\eta_i$. The coefficient of $\eta_i^2$ is $A=1/\tau^2+\sum_j1/\sigma_j^2$, and the linear term is $-2B\eta_i$, where $B=\mu_0/\tau^2+\sum_j z_j/\sigma_j^2$. Completing the square gives $A(\eta_i-B/A)^2$, plus a term that does not depend on $\eta_i$. Thus the normalized posterior is normal with variance $1/A$ and mean $B/A$:

\eta_i\mid z_1,\ldots,z_m\sim N\left( \frac{\mu_0/\tau^2+\sum_{j=1}^m z_j/\sigma_j^2} {1/\tau^2+\sum_{j=1}^m1/\sigma_j^2}, \frac{1}{1/\tau^2+\sum_{j=1}^m1/\sigma_j^2} \right).

This derivation depends on a normal prior and a linear Gaussian observation model in the parameter. It is not a general formula for nonlinear concentration models. In an NLME application, it can nevertheless provide an exact calculation when the chosen parameterization and observations have this form, and a useful check on numerical implementations.

Worked Example: One Gaussian Concentration Measurement

Suppose $\eta_i$ is an individual log-concentration level in a simplified model. Its population prior is $N(2.00,0.40^2)$, and one measured log concentration is $z=2.30$, with Gaussian residual standard deviation $0.20$. The model is $z=\eta_i+\epsilon$. Both the prior and residual model are normal, and the observation is linear in $\eta_i$, so the exact update above applies.

The prior precision is $1/0.40^2=6.25$, and the observation precision is $1/0.20^2=25.0$. Therefore the posterior variance is $1/(6.25+25.0)=0.0320$, and the posterior standard deviation is $\sqrt{0.0320}=0.179$.

The posterior mean is \[ \frac{(2.00)(6.25)+(2.30)(25.0)}{6.25+25.0} =\frac{12.5+57.5}{31.25}=2.24. \] Thus $\eta_i\mid z\sim N(2.24,0.0320)$. The data move the estimate from the prior mean of 2.00 toward the observed value 2.30, while uncertainty falls from a standard deviation of 0.40 to 0.179. If $\eta_i$ is the log of a concentration level in mg/L, $\exp(2.24)=9.40$ mg/L is the posterior median of that level, not its posterior mean.

Worked Example: Combining Two Measurements With Different Precision

Now suppose $\eta_i$ is a log-scale individual exposure parameter with prior $N(0,0.30^2)$. After subtracting known dose- and time-dependent offsets from two log concentrations, the measurements are $z_1=0.20$ and $z_2=0.50$. Their residual standard deviations are $0.20$ and $0.30$, respectively. Assume independent Gaussian errors and the linear model $z_j=\eta_i+\epsilon_j$. The measurement with lower residual variance should contribute more information.

The prior precision is $1/0.30^2=11.111$. The two observation precisions are $1/0.20^2=25.0$ and $1/0.30^2=11.111$. Their sum with the prior precision is $47.222$, so the posterior variance is $1/47.222=0.02118$, and the posterior standard deviation is $0.1455$.

The posterior mean is \[ \frac{(0)(11.111)+(0.20)(25.0)+(0.50)(11.111)}{47.222} =\frac{10.556}{47.222}=0.2235. \] Hence the posterior is $N(0.2235,0.02118)$. The result is closer to $0.20$ than to $0.50$, because the first measurement has greater precision. Both measurements contribute; simply averaging them would ignore their different error variances and the population prior.

Nonlinear Concentration Models

In a pharmacokinetic model, the expected concentration is often nonlinear in the individual parameters. For example, a one-compartment elimination model with fixed volume $V$ can have $f(t;k)=D/V\exp(-kt)$, where $k$ is an individual elimination-rate constant. With Gaussian additive concentration error, the posterior density is proportional to the population prior density times the Gaussian likelihood evaluated at this nonlinear prediction. There is generally no precision-weighted normal update.

Worked Example: Nonlinear Likelihood Updating on a Parameter Grid

For a transparent numerical illustration, use $C(t;k)=10\exp(-kt)$ mg/L, with time in hours and $k$ in inverse hours. At $t=2$ hours, the observed concentration is $3.00$ mg/L, and additive Gaussian residual error has known standard deviation $0.50$ mg/L. Represent the population prior on a three-point grid: $k=0.30,0.50,0.70$ h$^{-1}$, with prior probabilities $0.25,0.50,0.25$. This finite prior is a small-grid approximation for illustrating the likelihood weighting, not a claim that a continuous population distribution has only three possible values.

For each grid value, calculate the predicted concentration and the likelihood factor after dropping the Gaussian normalizing constant, which is common to all three values: $\exp\{-(3.00-C(2;k))^2/(2(0.50)^2)\}$. At $k=0.30$, $C=10\exp(-0.60)=5.488$ mg/L, so the factor is $\exp[-(3.00-5.488)^2/0.50]\approx0.0000042$. At $k=0.50$, $C=3.679$ mg/L, giving $\exp[-(3.00-3.679)^2/0.50]\approx0.398$. At $k=0.70$, $C=2.466$ mg/L, giving $\exp[-(3.00-2.466)^2/0.50]\approx0.565$.

Multiply each factor by its prior probability. The unnormalized posterior weights are approximately $0.0000011,\ 0.199,\ 0.141$; their sum is $0.340$. Dividing each by this sum gives posterior probabilities approximately $0.000003,\ 0.585,\ 0.415$, respectively. The observation strongly downweights the slowest elimination rate because it predicts a concentration substantially above the observed value. A continuous-prior analysis uses the same prior-times-likelihood principle but integrates over the parameter, typically with numerical integration or sampling.

Interpretation, Diagnostics, and Alternatives

The posterior expresses relative support for parameter values after combining the specified population distribution and the observed data model. It is conditional on those assumptions and on the values of $\psi$ used. A posterior mean, median, mode, or interval summarizes that distribution in a different way; none should be substituted for the distribution when uncertainty matters.

Several checks help identify a problematic update:

Check the observation model. Plot or inspect residuals against time and predicted concentration. A systematic pattern, changing spread, or correlation can indicate that the assumed independent constant-variance errors are inadequate.
Check prior-likelihood conflict. If the data favor parameter values in a region with negligible prior probability, the posterior may be pushed against a boundary or become highly sensitive to the prior. Verify the dose history, units, sampling times, covariates, and structural model before interpreting such a result.
Check posterior shape and computation. For nonlinear models, inspect the density or numerical draws for skewness, multiple modes, and boundary accumulation. A single normal approximation can hide these features; numerical integration or sampling must adequately explore the supported region.
Check what is treated as fixed. If $\psi$ and $\sigma$ are plugged in, the resulting individual posterior does not include uncertainty in those population and residual parameters.

A prior-only estimate ignores the individual concentrations. A likelihood-only estimate, such as an individual maximum-likelihood estimate, ignores the population distribution and can be unstable with sparse data. A posterior combines both sources and preserves uncertainty, while a point estimate such as a posterior mode compresses the distribution to one value. These approaches answer different questions; a point estimate is not a substitute for an individual posterior distribution.

The update is limited by the quality of the structural and residual models, the informativeness of the sampling schedule, and the population prior. Sparse or poorly timed samples may leave several parameter combinations plausible. A narrow posterior can also be misleading if the model is wrong or the prior is overconfident. For censored concentrations, nonlinear residual error, or correlated measurements, use an appropriate likelihood rather than forcing the Gaussian independent-error calculation.

Common pitfall: “The observed concentration is high, so the individual clearance must be low.” This reasoning treats one observation as a direct parameter measurement. A concentration depends on the full dose history, sampling time, other pharmacokinetic parameters, and residual error. The correction is to evaluate the likelihood under the complete concentration model and combine it with the population prior; the data update plausibility but need not identify one parameter uniquely.

Key takeaway: The individual posterior is proportional to the population prior times the likelihood of that individual’s observed concentrations. Precision-weighted normal formulas are exact only under their normal-linear assumptions; nonlinear pharmacokinetic models require the same Bayesian update with an appropriate numerical calculation.

Check Your Understanding

Answer each question, then compare with the solutions below.

If a prior is $N(1,0.5^2)$ and one linear Gaussian observation is $1.4$ with residual standard deviation $0.5$, what are the posterior mean and variance?
In the two-measurement example, why did the first concentration-derived value have more influence than the second?
For a nonlinear concentration model, is multiplying the prior density by the likelihood still valid if the posterior is not normal?
Name two assumptions that must be reconsidered if residuals show increasing variability with predicted concentration.

Solutions. (1) Both precisions are $1/0.5^2=4$, so the posterior variance is $1/(4+4)=0.125$ and the mean is $[(1)(4)+(1.4)(4)]/8=1.20$. (2) Its residual variance was $0.20^2$, smaller than $0.30^2$, so its precision was larger: 25 versus 11.111. (3) Yes. Bayes’ rule remains valid; the posterior simply may require numerical rather than closed-form calculation. (4) Reconsider constant residual variance and the choice of an additive Gaussian error model; a proportional or otherwise heteroscedastic error model may be more suitable.