Hierarchical Generative Models for Pharmacokinetic SDEs

Four Levels of a Mixed-Effects SDE

A pharmacokinetic SDE describes stochastic evolution within an individual. A mixed-effects SDE adds variation between individuals, then connects each latent trajectory to the measurements made on that individual. Keeping these levels distinct is essential: a population parameter is not a patient-specific random effect, a random effect is not a latent path, and neither is an observation.

Consider \(N\) individuals, with subject \(i\) observed at times \(t_{ij}\), for \(j=1,\ldots,n_i\). A useful generative description moves from population quantities to an individual effect, from that effect to an SDE path, and finally from the path to observations.

Definition: A hierarchical generative model specifies a sequence of conditional probability laws for the data. In a mixed-effects SDE, these levels are: (1) population parameters \(\theta\), which govern the population distributions and model functions; (2) individual random effects \(\eta_i\), drawn from a distribution indexed by population parameters; (3) the latent process \(X_i(t)\), generated by an SDE conditional on the individual effect and known inputs; and (4) observations \(Y_{ij}\), generated conditional on the latent state at each measurement time and any relevant individual quantities.

A common parameterization maps the individual effect to subject-specific parameters by \(\phi_i=h(\beta,\eta_i)\). Here \(\beta\) is a population-level parameter and \(h\) is a specified link or transformation. For example, \(CL_i=CL_{\mathrm{pop}}\exp(\eta_i)\) makes clearance positive. The random effect \(\eta_i\) is a random variable in the population model; after an individual has been sampled, its realized value is a fixed but usually unknown subject-specific quantity.

Conditional on \(\phi_i\) and known inputs \(u_i(t)\), the latent state follows a model such as

dX_i(t)=b\bigl(X_i(t),t,u_i(t);\phi_i\bigr)\,dt +G\bigl(X_i(t),t,u_i(t);\phi_i\bigr)\,dW_i(t).

The observation law might be written \(Y_{ij}\mid X_i(t_{ij}),\phi_i\sim g(\,\cdot\mid X_i(t_{ij}),\phi_i)\). It may represent assay error, censoring, or another measurement process. The process noise is the randomness generated by the SDE; observation error is generated by the observation law. They are different sources of variation, even when both affect observed concentrations.

Assumptions That Make the Hierarchy Well Defined

A model needs conditions that make its stated probability laws meaningful. One standard set of sufficient conditions is that, for each allowed \(\phi_i\) and known input, \(b\) and \(G\) are measurable in time, globally Lipschitz in the state, and satisfy a linear-growth bound. These conditions give a unique nonexplosive strong solution for a specified initial state with finite second moment. More general existence conditions are possible, but should be stated when used.

The random effects have a specified population distribution \(F_\eta(\,\cdot\mid\theta)\). Across individuals they are usually assumed independent and identically distributed. The subject-specific Brownian motions are independent across individuals and independent of the random effects, conditional on population parameters. Within an individual, observations are commonly assumed conditionally independent given the latent path and individual quantities. This last assumption does not mean the observations are independent after the path and effects are integrated out.

The population parameter \(\theta\) is treated as fixed in the sampling model. It may be unknown and estimated from data. If a Bayesian analysis assigns it a prior, that introduces an additional probability level; it does not change the meaning of the conditional generative model given \(\theta\).

Worked Example: Assigning a Clearance Effect to the Right Level

Suppose a one-compartment amount \(A_i(t)\) has known volume \(V\), known dose history \(u_i(t)\), and individual clearance \(CL_i\). Let

CL_i=CL_{\mathrm{pop}}\exp(\eta_i),\qquad dA_i(t)=\left[-\frac{CL_i}{V}A_i(t)+u_i(t)\right]dt+\sigma\,dW_i(t), \qquad Y_{ij}=\frac{A_i(t_{ij})}{V}+\epsilon_{ij}.

Take \(CL_{\mathrm{pop}}=5\) L/h, \(V=10\) L, and \(\eta_i\sim N(0,0.04)\). These are population-level specifications: \(CL_{\mathrm{pop}}\), \(V\), \(\sigma\), and the variance \(0.04\) are parameters. For a realized effect \(\eta_i=0.1\), the individual clearance is \(CL_i=5e^{0.1}\approx5.526\) L/h, so the individual elimination rate is \(CL_i/V\approx0.553\) h\(^{-1}\).

The amount \(A_i(t)\) is the latent process, driven by its SDE and the known dose input. The observed concentration \(Y_{ij}\) is a separate random variable; for example, one could specify \(\epsilon_{ij}\sim N(0,0.01)\), independently across measurements conditional on the path. Thus \(\eta_i\) explains between-subject variation in clearance, \(\sigma\) scales within-subject process noise, and the variance of \(\epsilon_{ij}\) scales measurement error. They cannot be substituted for one another without changing the model.

From the Generative Model to the Data Distribution

For a subject, let \(y_i=(y_{i1},\ldots,y_{in_i})\), and let \(P_{\phi_i}(dx)\) denote the probability law of the latent path conditional on \(\phi_i\) and the known inputs. Suppose the observation density \(g(y\mid x,\phi)\) is normalized and measurable, and that the conditional observation errors are independent given the path and individual quantities. The subject’s marginal density is then

p(y_i\mid\theta) = \int \left[ \int \prod_{j=1}^{n_i} g\bigl(y_{ij}\mid x(t_{ij}),\phi\bigr) \,P_{\phi}(dx) \right] F_{\eta}(d\eta\mid\theta), \qquad \phi=h(\beta,\eta).

This expression integrates over two distinct unknown sources: the random effect and the latent SDE path. It makes no claim that either integral has a closed form. In particular, replacing the random effect by its mean or the path by its conditional mean generally changes the data distribution.

Result: Under the stated random-effect and path independence assumptions across subjects, and conditional on fixed \(\theta\), the joint observed-data density factorizes as \(p(y_1,\ldots,y_N\mid\theta)=\prod_{i=1}^{N}p(y_i\mid\theta)\), with each factor given by the marginal expression above.

Condition on the individual effect.
For subject \(i\), fix \(\eta_i=\eta\), and hence \(\phi_i=h(\beta,\eta)\). The SDE then defines the conditional path law \(P_\phi\).

Condition on the path.
By conditional independence of observations, the conditional density of \(y_i\) is \(\prod_j g(y_{ij}\mid x(t_{ij}),\phi)\).

Average over the path and effect.
The law of total probability first integrates this product with respect to \(P_\phi(dx)\), then integrates with respect to \(F_\eta(d\eta\mid\theta)\), giving \(p(y_i\mid\theta)\).

Use independence across subjects.
The random effects, Brownian drivers, and observation mechanisms are independent across subjects under the model. Thus the full joint density is the product of the subject-level marginal densities. \(\square\)

This factorization is the basis for a likelihood built from independent subjects. It does not permit factorization over measurements within one subject: those measurements share both an individual effect and a latent path.

Worked Example: A Conditional OU Trajectory and Its Measurement

Let \(Z_i(t)\) be a latent log-concentration state satisfying

dZ_i(t)=-k_iZ_i(t)\,dt+0.1\,dW_i(t),\qquad Z_i(0)=1,\qquad k_i=0.5e^{\eta_i},

with \(\eta_i\sim N(0,0.04)\), and let \(Y_i(1)=Z_i(1)+\epsilon_i\), where \(\epsilon_i\sim N(0,0.05)\) independently of the path and effect. Conditional on \(\eta_i\), the Ornstein–Uhlenbeck solution has

E[Z_i(1)\mid\eta_i]=e^{-k_i}, \qquad \operatorname{Var}(Z_i(1)\mid\eta_i) =\frac{0.1^2}{2k_i}\left(1-e^{-2k_i}\right).

For a realized value \(\eta_i=0.2\), \(k_i=0.5e^{0.2}\approx0.6107\). Therefore \(E[Z_i(1)\mid\eta_i]\approx e^{-0.6107}\approx0.543\), and \[ \operatorname{Var}(Z_i(1)\mid\eta_i) \approx \frac{0.01}{1.2214}(1-e^{-1.2214}) \approx \frac{0.01}{1.2214}(1-0.2949) \approx0.00577. \] The conditional observation variance is \(0.00577+0.05=0.05577\), since the independent measurement-error variance adds to the latent-state variance. These are conditional calculations for the specified effect, not population-marginal moments; obtaining those requires averaging over \(\eta_i\) as well.

Why Measurements From One Individual Are Dependent

Even conditionally independent observation errors do not make observations independent after the latent quantities are marginalized. Let \(H_i\) denote the subject’s random effect and entire latent path. Suppose \(Y_{i1}\) and \(Y_{i2}\) are conditionally independent given \(H_i\), have finite second moments, and write \(m_j(H_i)=E[Y_{ij}\mid H_i]\). The law of total covariance gives

\operatorname{Cov}(Y_{i1},Y_{i2}) = E[\operatorname{Cov}(Y_{i1},Y_{i2}\mid H_i)] + \operatorname{Cov}(m_1(H_i),m_2(H_i)) = \operatorname{Cov}(m_1(H_i),m_2(H_i)).

For completeness, expand the covariance using conditional expectation: \(E[Y_{i1}Y_{i2}]=E[E[Y_{i1}Y_{i2}\mid H_i]]\). Conditional independence makes the inner expectation \(m_1(H_i)m_2(H_i)\). Also \(E[Y_{ij}]=E[m_j(H_i)]\). Subtracting the product of these means yields precisely \(\operatorname{Cov}(m_1(H_i),m_2(H_i))\). Thus a shared random effect or correlated latent trajectory can induce marginal dependence even when the measurement errors are independent. \(\square\)

Worked Example: Shared Random Effect Creates Correlated Measurements

Consider the simple latent SDE \(dZ_i(t)=(\beta+\eta_i)dt+0.2\,dW_i(t)\), \(Z_i(0)=0\), with \(\eta_i\sim N(0,0.25)\). Observe \(Y_i(t)=Z_i(t)+\epsilon_i(t)\), with independent errors of variance \(0.01\). At times 1 and 2, the random-effect contribution to each latent state is \(t\eta_i\), and the Brownian contribution is \(0.2W_i(t)\).

The shared-effect covariance is \(1\cdot2\cdot0.25=0.50\). Since \(\operatorname{Cov}(W_i(1),W_i(2))=1\), the shared Brownian contribution is \(0.2^2\cdot1=0.04\). The independent measurement errors contribute no cross-time covariance. Hence \(\operatorname{Cov}(Y_i(1),Y_i(2))=0.50+0.04=0.54\).

The variances are \(\operatorname{Var}(Y_i(1))=1^2(0.25)+0.2^2(1)+0.01=0.30\), and \(\operatorname{Var}(Y_i(2))=2^2(0.25)+0.2^2(2)+0.01=1.09\). Thus the correlation is \(0.54/\sqrt{0.30\cdot1.09}\approx0.944\). Treating the two measurements as marginally independent would discard substantial dependence generated by the hierarchy.

Model Checks, Alternatives, and Limits

A practical specification check is to label every quantity by its level before fitting: fixed population parameter, random effect, latent state or path, known input, or observation. Then check the conditional laws in the intended order. Simulated subjects should display between-subject variation consistent with the random-effect distribution, within-subject variation consistent with the SDE, and measurement scatter consistent with the observation model. If simulated trajectories vary when only measurement error is supposed to vary, or if independent subjects share a Brownian driver unintentionally, the implementation does not match the stated hierarchy.

A mixed-effects SDE differs from two common alternatives. A population-average SDE with no random effects describes within-process stochasticity but does not represent persistent between-individual parameter variation. A subject-specific fixed-parameter analysis assigns each individual a separate unknown parameter rather than drawing effects from a population distribution; it may be useful for individual inference but does not by itself specify a population distribution for new individuals. A deterministic ODE with observation error alone omits process noise. These models answer different questions and are not interchangeable merely because they can produce similar fitted means.

A hierarchical specification also does not guarantee that its parameters can be estimated from the available data. Different combinations of random-effect variation, process noise, and observation error may produce similar observed variability. Sparse sampling can make these components difficult to distinguish. The identifiability material and the earlier material on Stochastic Calculus provide relevant foundations, but the hierarchy must still be checked against the actual observation schedule and measurement design. A model can be mathematically well defined yet weakly informed by data.

Common pitfall: It is incorrect to say, “The observations are independent because the assay errors are independent.” That reasoning conditions on neither the shared random effect nor the latent path. The correction is to state the conditioning explicitly: observations may be independent given the effect and path, while being dependent under the marginal data distribution.

Check Your Understanding

Answer each question by identifying the relevant level of the hierarchy. Brief solutions follow.

In \(CL_i=CL_{\mathrm{pop}}\exp(\eta_i)\), which quantities are population-level and which vary by individual?
Why does an SDE diffusion coefficient not generally represent assay error?
Under what independence assumptions does the full observed-data density factorize across subjects?
Can conditionally independent measurements from one subject be marginally dependent? Give the reason.
What must be integrated out to obtain a subject’s marginal data density?

Solutions. (1) \(CL_{\mathrm{pop}}\) is a population parameter; \(\eta_i\) is an individual random effect, and \(CL_i\) is the resulting individual parameter. (2) The diffusion coefficient governs randomness in the latent path; assay error enters through the observation law. (3) Subject effects, Brownian drivers, and observation mechanisms must be independent across subjects conditional on fixed population parameters. (4) Yes. A shared effect or latent path can make conditional means covary after marginalization. (5) Integrate over both the random effect and the latent path, using their specified conditional laws.