Hierarchical Models For Formulation Comparisons

Why Use a Hierarchical Model?

Concentration-time profiles from different people rarely follow exactly the same curve. A person may have higher clearance, a different volume of distribution, or a different absorption rate than another person. At the same time, two formulations may systematically produce different concentration profiles. A useful model must represent both patterns without treating every observed concentration as an independent measure of the formulation effect.

A nonlinear mixed-effects model does this by combining a population model, subject-level random effects, and a residual-error model. The population model describes typical PK parameters and the average formulation contrast. Random effects describe how a subject’s parameters differ from those population values. Residual error describes the remaining observation-level variation after accounting for the model and subject.

Definition: A hierarchical nonlinear mixed-effects model specifies a nonlinear mean function for each observation, lets one or more model parameters vary across subjects through subject-level random effects, and describes residual variation conditional on those subject effects. Population parameters describe the distribution or typical values of the subject-specific parameters; they are not themselves individual measurements.

Model Structure for a Formulation Comparison

Suppose subject \(i\) contributes concentration measurements under reference (R) and test (T) formulations. Let \(t_{ijk}\) be the time of measurement \(k\) for subject \(i\) under formulation \(j\), and let \(z_{ij}=0\) for R and \(z_{ij}=1\) for T. The formulation contrast can enter a PK parameter as a fixed population effect. For example, use a one-compartment oral model with first-order absorption, and let the formulation multiply the subject’s baseline bioavailability scale:

F_{ij}=F_i\exp(\beta z_{ij}), \qquad C_{ij}^{\mathrm{pred}}(t)= \frac{D_{ij}F_{ij}k_{a,i}}{V_i(k_{a,i}-k_{e,i})} \left[\exp(-k_{e,i}t)-\exp(-k_{a,i}t)\right], \qquad k_{e,i}=\frac{CL_i}{V_i}.

Here \(D_{ij}\) is the administered dose, \(CL_i\) is clearance, \(V_i\) is volume, and \(k_{a,i}\) is the absorption-rate constant. The population formulation contrast is \(\beta\): the test-to-reference multiplier on this model scale is \(\exp(\beta)\). The scale \(F_i\) is the subject’s reference-formulation scale. In oral data, an absolute bioavailability parameter may not be separately identifiable from clearance; the formulation multiplier should therefore be interpreted according to the model and design, not automatically as an absolute bioavailability measurement.

A common population model for the subject parameters is:

\begin{aligned} \log(CL_i)&=\log(CL_{\mathrm{pop}})+\eta_{CL,i},\\ \log(V_i)&=\log(V_{\mathrm{pop}})+\eta_{V,i},\\ \log(k_{a,i})&=\log(k_{a,\mathrm{pop}})+\eta_{ka,i},\\ \log(F_i)&=\log(F_{\mathrm{pop}})+\eta_{F,i}. \end{aligned}

The vector of random effects for subject \(i\), denoted \(\boldsymbol{\eta}_i\), is commonly modeled as multivariate normal with mean zero and covariance matrix \(\Omega\). A nonzero off-diagonal element of \(\Omega\) represents covariance between two subject-level effects. Exponentiating a normally distributed random effect gives positive individual parameters. These are assumptions to assess, not facts guaranteed by using this model.

In this formulation, \(\beta\) is a population-level fixed contrast, while the \(\eta\) terms are subject-level deviations. The same subject effect is shared across that subject’s observations and, in a crossover, across their formulation periods. This model does not add a separate random subject-by-formulation effect; doing so would represent a different source of variation and requires its own design and interpretation.

Worked Example: Specifying the Formulation Contrast

Consider a hypothetical study in which the estimated population contrast is \(\beta=0.1823\). The model places the formulation effect on the log scale, so the test-to-reference multiplier is:

\exp(\beta)=\exp(0.1823)=1.200.

Thus, under this model, the test formulation multiplies the reference formulation’s bioavailability-related scale by 1.20 at the population level, or 20.0% higher. This interpretation concerns the parameter where the contrast was placed. It is not, by itself, a statement that every subject has a 20.0% higher value or that a derived exposure metric has exactly the same ratio.

For illustration, suppose the reference scale for a hypothetical subject is \(F_i=0.500\). The model gives \(F_{iT}=0.500(1.200)=0.600\). This calculation illustrates the multiplicative contrast; it does not establish that an unconstrained model parameter is a physically measured absolute fraction absorbed.

Interindividual Variability and Residual Error

For positive concentrations, one possible residual model is lognormal:

\log(C_{ijk}^{\mathrm{obs}})= \log\!\left(C_{ij}^{\mathrm{pred}}(t_{ijk})\right)+\epsilon_{ijk}, \qquad \epsilon_{ijk}\sim N(0,\sigma^2).

The residual errors are often assumed conditionally independent given the subject effects and model covariates, with mean zero and common variance \(\sigma^2\). Under this specification, \(\sigma\) describes observation-level error on the log-concentration scale. A residual standard deviation of 0.10 is not an error of 0.10 concentration units; it is a log-scale standard deviation. Other residual structures, such as combined additive and proportional error, may be more suitable when the assay and concentration range support them.

The three levels have distinct meanings. Population parameters such as \(CL_{\mathrm{pop}}\) and \(\beta\) describe the typical model and average formulation contrast. Random effects such as \(\eta_{CL,i}\) describe persistent between-subject differences in parameters. Residual errors \(\epsilon_{ijk}\) describe observation-level deviations remaining after the predicted curve and subject effects are accounted for. Repeated measurements from the same person are linked through the shared random effects, so they are not equivalent to independent subjects.

Interpretation: In a lognormal model, a fixed effect is a population-level shift on the log scale; exponentiating it gives a multiplicative contrast. The covariance matrix \(\Omega\) describes interindividual variability and association among subject-level parameters. The residual variance \(\sigma^2\) describes conditional observation-level scatter. These quantities are not interchangeable.

Why the Hierarchy Matters: A Short Derivation

The variance decomposition shows why ignoring subject effects can misrepresent the data. Consider a simplified log-scale observation \(Y_{ik}=m+\eta_i+\epsilon_{ik}\), where \(m\) is a population mean, \(\eta_i\) is a subject effect with variance \(\omega^2\), and \(\epsilon_{ik}\) is residual error with variance \(\sigma^2\). Assume the effects are independent and both have mean zero. Then:

\begin{aligned} \operatorname{Var}(Y_{ik}) &=\operatorname{Var}(\eta_i+\epsilon_{ik})\\ &=\operatorname{Var}(\eta_i)+\operatorname{Var}(\epsilon_{ik})\\ &=\omega^2+\sigma^2. \end{aligned}

For two observations from the same subject, with independent residual errors, the covariance is instead:

\operatorname{Cov}(Y_{ik},Y_{i\ell}) =\operatorname{Cov}(\eta_i+\epsilon_{ik},\eta_i+\epsilon_{i\ell}) =\omega^2,\qquad k\ne\ell.

The shared subject effect creates within-subject dependence. A model that assigns all variation to independent residual errors misses this covariance structure. In a nonlinear PK model, the exact marginal distribution is more involved because the concentration function is nonlinear in the random effects, but the same hierarchy still determines which observations share a subject-level source of variation.

Worked Example: Translating Random Effects Into Individual Parameters

Suppose a hypothetical population model has \(CL_{\mathrm{pop}}=5.00\) L/h and \(V_{\mathrm{pop}}=50.0\) L. For one hypothetical subject, let \(\eta_{CL,i}=0.100\) and \(\eta_{V,i}=-0.0500\). The individual parameters are:

CL_i=5.00\exp(0.100)=5.526\ \text{L/h}, \qquad V_i=50.0\exp(-0.0500)=47.56\ \text{L}.

The individual elimination rate is therefore:

k_{e,i}=\frac{5.526\ \text{L/h}}{47.56\ \text{L}} =0.1162\ \text{h}^{-1}.

These are individual parameter values obtained by combining population estimates with this subject’s random effects. They do not mean that the population parameters themselves changed. If the random effects for log clearance and log volume have standard deviations 0.200 and 0.150 and correlation 0.500, their covariance is \(0.500(0.200)(0.150)=0.0150\). The corresponding covariance matrix is positive definite because its determinant is \(0.0400(0.0225)-(0.0150)^2=0.000675>0\), a necessary consistency check for this two-effect covariance matrix.

Worked Example: Checking a Conditional Residual

Suppose the model predicts a concentration of \(0.7593\) mg/L at a particular time for a hypothetical observation, and the observed concentration is \(0.8000\) mg/L. Under the lognormal residual model, the residual is:

r=\log(0.8000)-\log(0.7593) =\log\left(\frac{0.8000}{0.7593}\right) \approx 0.0522.

If the fitted residual standard deviation is \(\sigma=0.100\), the standardized conditional residual is \(r/\sigma=0.522\). This observation is about 0.522 residual standard deviations above its prediction on the log scale. It is not a deviation of 0.522 mg/L.

For a simple log-scale model with a random-effect standard deviation of 0.120 and residual standard deviation of 0.100, the marginal log-scale variance is \(0.120^2+0.100^2=0.0244\), and its standard deviation is \(\sqrt{0.0244}=0.1562\). This calculation separates the two sources of variation; it does not imply that the nonlinear concentration data have a constant-scale normal distribution.

Assumptions, Diagnostics, and Alternatives

A fitted hierarchical model is interpretable only to the extent that its structural and statistical assumptions are adequate. Important checks include whether the concentration-time function captures the observed curve shape, whether the random-effects distribution is plausible, whether residual spread changes with concentration or time, and whether the model can identify the parameters being estimated. The earlier material on identifiability of bioavailability and clearance is especially relevant: adding random effects does not resolve a structural confounding that the data cannot distinguish.

Structural mismatch: systematic residual patterns across time, such as underprediction during absorption and overprediction later, can indicate a poorly chosen structural model or an omitted process.
Residual misspecification: residual spread that grows with concentration, time trends, or strong outliers can contradict a constant log-scale residual model. Consider whether a different residual structure is scientifically and statistically justified.
Random-effects problems: implausible variance estimates, correlations near their boundaries, or unstable estimates can indicate weak information, over-parameterization, or an unsuitable random-effects structure.
Identifiability and design: sparse sampling, too few subjects, or limited observations around informative portions of the curve can make parameters difficult to distinguish. A converged optimizer does not prove that parameters are well identified.
Influential subjects: inspect subject-level predictions and residuals. A small number of unusual profiles can dominate estimates of variability or the formulation contrast.

A fixed-effects-only nonlinear model is simpler but assumes that subjects share the modeled parameters apart from residual error; it can understate dependence and fail to represent real interindividual variability. Fitting separate curves to each subject avoids a common population distribution but can be unstable with sparse individual data and does not efficiently estimate a shared population contrast. A hierarchical model partially pools information: it estimates population parameters while allowing subjects to differ. Its additional complexity is useful only when supported by the design and data.

The model’s limitations should remain explicit. Distributional assumptions for random effects may be imperfect, and estimates can be sensitive to the structural model, residual specification, and sampling schedule. A population formulation contrast is an average model parameter, not a claim that every individual responds identically. Nor does the model alone establish a regulatory bioequivalence conclusion: the relevant estimand, exposure measure, uncertainty calculation, and decision rule must also be specified.

Common Mistake

Pitfall: “The formulation effect is \(\beta\), so \(\beta\) is the percent increase in concentration for every subject.” This is incorrect. In the example model, \(\beta\) is a population log-scale contrast on the parameter where the formulation effect was placed. Its multiplicative interpretation is \(\exp(\beta)\), while individual profiles also depend on subject parameters, time, and residual variation. Report the contrast on its stated scale and do not equate it automatically with every subject’s concentration or a derived exposure ratio.

Check Your Understanding

Answer each question, then compare your reasoning with the short solutions below.

If a fitted formulation coefficient is \(\beta=0.0953\), what is its multiplicative interpretation on the modeled parameter scale?
In the stated model, why do two concentrations from the same subject generally have correlated deviations from the population curve?
What does \(\sigma=0.08\) represent in a lognormal residual model, and what does it not represent?
A model estimates a random-effects standard deviation of zero for one parameter. Does that prove that no subjects differ in that parameter?
Name one residual diagnostic pattern that could suggest model misspecification and describe what it may indicate.

Solutions. (1) \(\exp(0.0953)=1.100\), so the modeled parameter is multiplied by approximately 1.10, or 10.0% higher, for T than R at the population level. (2) They share the same subject-level random effects, which induce dependence across that subject’s observations. (3) It is a residual standard deviation on the log-concentration scale; it is not an error of 0.08 concentration units. (4) No. A boundary estimate can reflect limited information, model structure, or estimation constraints; it is not proof of identical subjects. (5) For example, residuals that are systematically positive early and negative late may indicate that the structural concentration-time model does not capture the observed time course.