End-to-End Estimation of a Mixed-Effects SDE

From a population model to an estimable likelihood

A mixed-effects SDE describes two sources of subject-to-subject variation at once: differences in typical pharmacokinetic parameters and stochastic variation along each subject’s concentration trajectory. Estimation must account for both. In particular, a subject’s latent random effect is shared across that subject’s observations, whereas the Brownian increments and measurement errors vary over time. Treating all these terms as independent observation-level noise gives the wrong likelihood and can distort estimates of both variability and elimination.

We will work with a one-compartment elimination model on the log-concentration scale. Let $C_i(t)>0$ be the concentration in subject $i$, and let $X_i(t)=\log C_i(t)$. A normally distributed random effect changes the subject’s elimination rate. Observations are log concentrations with additive measurement error. This setup keeps the concentration positive, makes the role of each variance component explicit, and permits a transparent likelihood calculation.

Definition: For subjects $i=1,\ldots,n$, consider $$ dC_i(t)=-k\exp(\eta_i)C_i(t)\,dt+\sigma C_i(t)\,dW_i(t),\qquad C_i(0)=C_0, $$ with observations $$ Z_{ij}=X_i(t_{ij})+\epsilon_{ij},\qquad X_i(t)=\log C_i(t). $$ Here $k>0$ is the population elimination-rate parameter, $\sigma>0$ is the diffusion scale, $\eta_i\sim N(0,\omega^2)$ is a subject-specific random effect, and $\epsilon_{ij}\sim N(0,\tau^2)$ is measurement error. The Brownian motions, random effects, and errors are mutually independent across subjects and sources; observation times are treated as fixed. The parameters to estimate are $\theta=(k,\sigma,\omega,\tau)$, with positive components.

This model presumes that a single elimination rate is meaningful over the observation window, that the random effects are drawn from the stated population distribution, and that measurement errors are independent with constant variance on the log scale. These are modelling assumptions, not consequences of the SDE. The earlier material on hierarchical generative models and nested integration supplies the broader distinction between population variables and latent paths; here the goal is to turn that specification into an estimation procedure.

Conditional trajectory and observation likelihood

Apply Itô’s formula to $X_i(t)=\log C_i(t)$. If $dC_i=aC_i\,dt+bC_i\,dW_i$, then $d\log C_i=(a-b^2/2)dt+b\,dW_i$. In this model $a=-k\exp(\eta_i)$ and $b=\sigma$, so

dX_i(t)=-\left(k\exp(\eta_i)+\frac{\sigma^2}{2}\right)dt+\sigma\,dW_i(t).

For fixed $\eta_i$, integration from time $0$ to $t$ gives a Gaussian process with mean $X_i(0)-(k\exp(\eta_i)+\sigma^2/2)t$ and Brownian covariance $\sigma^2\min(s,t)$. Adding independent observation error therefore gives a multivariate normal observation vector conditional on $\eta_i$.

Conditional observation distribution: Write $Z_i=(Z_{i1},\ldots,Z_{im_i})^\mathsf{T}$. Conditional on $\eta_i=\eta$, $$ Z_i\mid\eta\sim N\!\left(\mu_i(\eta),\Sigma_i\right), $$ where $$ \mu_{ij}(\eta)=\log C_0-\left(k\exp(\eta)+\frac{\sigma^2}{2}\right)t_{ij}, \qquad (\Sigma_i)_{j\ell}=\sigma^2\min(t_{ij},t_{i\ell})+\tau^2\mathbf{1}_{\{j=\ell\}}. $$

The shared random effect enters the conditional mean at every time. The Brownian term instead contributes covariance between two observations in proportion to the time they share before the earlier observation; measurement error contributes only to the diagonal. These distinctions are essential when constructing the likelihood.

Result: Subject $i$’s marginal likelihood contribution is $$ L_i(\theta)=\int_{-\infty}^{\infty} \phi\!\left(Z_i;\mu_i(\eta),\Sigma_i\right) \phi(\eta;0,\omega^2)\,d\eta, \qquad L(\theta)=\prod_{i=1}^{n}L_i(\theta), $$ where $\phi(\,\cdot\,;m,V)$ denotes a normal density with mean $m$ and covariance $V$. The integral is over the subject-level random effect; the latent Brownian path has already been integrated out in the conditional Gaussian distribution.

This expression follows by conditioning on $\eta_i$ and applying the law of total probability: the joint density of $Z_i$ is the integral of its conditional density times the density of $\eta_i$. Independence between subjects then gives the product. Thus the marginal log-likelihood is $\ell(\theta)=\sum_i\log L_i(\theta)$. In this model the conditional distribution is available exactly; in less tractable SDEs, transition densities or latent-path integrals may also need numerical approximation.

Variance components and likelihood approximation

For numerical optimization, estimate unconstrained parameters such as $\alpha_k=\log k$, $\alpha_\sigma=\log\sigma$, $\alpha_\omega=\log\omega$, and $\alpha_\tau=\log\tau$. Transforming back by exponentiation guarantees positive rate and scale estimates. For a vector random effect with covariance matrix $\Omega$, a standard extension is $\Omega=LL^\mathsf{T}$, with positive diagonal entries of the lower-triangular matrix $L$ parameterized on the log scale. This avoids proposing an invalid covariance matrix during optimization.

The integral over $\eta$ can be approximated by Gauss–Hermite quadrature. If $u_q,w_q$ are nodes and weights for integration against the standard normal density, then

L_i(\theta)\approx \sum_{q=1}^{Q}w_q\, \phi\!\left(Z_i;\mu_i(\omega u_q),\Sigma_i\right), \qquad \sum_{q=1}^{Q}w_q=1.

The nodes are scaled by $\omega$ because $\eta_i$ has standard deviation $\omega$. This is a deterministic approximation to the random-effect integral, not a replacement of each subject’s random effect by zero. Increasing $Q$, checking convergence of the objective and estimates, or comparing with adaptive quadrature are useful ways to detect inadequate integration.

Worked Example: A Three-Node Marginal Likelihood

Consider one subject observed at $t=1$, with $Z-\log C_0=-0.55$. Let $k=0.4$, $\sigma=0.2$, $\omega=0.3$, and $\tau=0.1$. At this single time the conditional variance is $\sigma^2+\tau^2=0.05$. Use the three-node standard-normal rule with nodes $u=0,\pm\sqrt{3}$ and weights $2/3,1/6,1/6$.

The corresponding random-effect values are $0$ and $\pm0.3\sqrt{3}=\pm0.5196$. The conditional mean of $Z-\log C_0$ is $-[0.4\exp(\eta)+0.02]$. Thus the means at $\eta=0,0.5196,-0.5196$ are approximately $-0.4200,-0.6927,-0.2579$, respectively. The observed-minus-mean residuals are $-0.1300,0.1427,-0.2921$.

A normal density with variance $0.05$ has prefactor $1/\sqrt{2\pi(0.05)}\approx1.7841$. The three conditional densities are therefore $$ 1.7841e^{-0.1300^2/0.1}\approx1.507,\qquad 1.7841e^{-0.1427^2/0.1}\approx1.455,\qquad 1.7841e^{-0.2921^2/0.1}\approx0.760. $$ The quadrature approximation to the marginal density is $$ \frac{2}{3}(1.507)+\frac{1}{6}(1.455)+\frac{1}{6}(0.760) \approx1.374. $$ This is a likelihood contribution for the observed log concentration under the stated parameter values. It is not yet a parameter estimate: estimation requires maximizing the sum of log contributions across all subjects.

Specify the hierarchy.
Choose the PK drift and diffusion, random-effect distribution, observation model, and independence assumptions before constructing the likelihood.

Form subject contributions.
Condition on each random effect, compute or approximate the observation density, and integrate over that effect’s population distribution.

Optimize and verify.
Maximize the marginal log-likelihood using positive-scale parameterizations, then check numerical integration accuracy and whether estimates are supported by the sampling design.

Estimation, variance components, and interpretation

A practical fit evaluates $\ell(\theta)$ at trial parameter values, uses a numerical optimizer on unconstrained transformed parameters, and repeats the calculation with a more accurate integration rule. Multiple initial values are prudent because likelihood surfaces can be flat or have local optima. The reported estimate $\hat{\theta}$ is the maximizer found; it should not be treated as reliable merely because an optimizer returned a convergence code.

The diffusion and assay-error components may be difficult to separate if observations are sparse or nearly simultaneous. The model’s covariance structure explains why: Brownian variation accumulates with time and induces correlation across observations, while independent measurement error increases only marginal variances. Separating a random effect from these sources also requires enough repeated observations per subject and enough subjects to characterize the population distribution.

Worked Example: Starting Values from a Covariance Decomposition

A useful starting-value calculation can be made for a simpler local model with a random intercept $\eta_i\sim N(0,\omega^2)$, Brownian variation $\sigma W_i(t)$, and independent measurement error of variance $\tau^2$. This is a covariance-based approximation for initialization, not the nonlinear random-clearance likelihood above. Suppose observations are at $t_1=1$ and $t_2=2$, and the empirical covariance matrix of centered observations is $$ S=\begin{pmatrix}0.14&0.11\\0.11&0.17\end{pmatrix}. $$

The model implies $\operatorname{Var}(Z(t))=\omega^2+\sigma^2t+\tau^2$ and $\operatorname{Cov}(Z(t_1),Z(t_2))=\omega^2+\sigma^2t_1$, since the shared random intercept and Brownian path contribute covariance, but independent errors do not. Subtracting the two variances gives $0.17-0.14=\sigma^2(2-1)$, so $\sigma^2=0.03$. The off-diagonal covariance gives $0.11=\omega^2+0.03$, hence $\omega^2=0.08$. Finally, $0.14=0.08+0.03+\tau^2$, so $\tau^2=0.03$. Substitution predicts the second variance $0.08+2(0.03)+0.03=0.17$, matching the supplied matrix.

These are moment-based starting values, not maximum-likelihood estimates, and empirical covariance entries can be noisy. A non-positive component estimate or a covariance pattern inconsistent with the model is a warning that the simple decomposition is inadequate or weakly informed by the data; it is not a reason to silently force a scientifically meaningful estimate.

Worked Example: Approximate Individual Random-Effect Prediction

Continue the single-observation calculation above. The three unnormalized quadrature contributions to the posterior weights are the prior quadrature weights times the conditional densities: for $\eta=0.5196,0,-0.5196$, these are $1.455/6=0.2425$, $2(1.507)/3=1.0047$, and $0.760/6=0.1267$. Their sum is approximately $1.374$.

Normalizing gives weights $0.1765,0.7312,0.0922$, which sum to approximately one. The quadrature approximation to the posterior mean of the random effect is therefore $$ E[\eta\mid Z]\approx(0.5196)(0.1765)+(0)(0.7312)+(-0.5196)(0.0922)\approx0.044. $$ The corresponding posterior mean elimination rate is approximately $$ E[k\exp(\eta)\mid Z]\approx0.4\{(1.6818)(0.1765)+(1)(0.7312)+(0.5948)(0.0922)\}\approx0.433. $$ This is an individual estimate conditional on the population parameters and this subject’s observation. It is not a new population estimate, and it should not be interpreted as if the posterior uncertainty had vanished.

Checks, alternatives, and limitations

A sound end-to-end analysis checks both the scientific model and the numerical calculation. Verify that the fitted diffusion and error scales are positive, the optimizer has not merely stopped at a boundary, and estimates are stable to starting values and quadrature accuracy. Examine whether the observed time grid can distinguish the proposed variance components. A likelihood that changes materially when the number of quadrature nodes increases is not yet numerically trustworthy.

Do not set every random effect to zero. Evaluating only the conditional density at $\eta=0$ ignores population heterogeneity; it is generally not the marginal likelihood.
Do not multiply conditional subject likelihoods before integrating a shared effect. For a subject, the effect is common to all observations, so integrate the joint density over one $\eta_i$, not separate effects for each time point.
Do not confuse Brownian variation with measurement error. Brownian increments create time-dependent cross-observation covariance; independent assay error does not.
Do not rely on a single integration setting or optimizer start. Sensitivity to either indicates that numerical or inferential uncertainty needs attention.

For the present model, Gaussian transition calculations integrate the latent path exactly, leaving only a one-dimensional random-effect integral. Gauss–Hermite quadrature is often effective for such integrals. Alternatives include adaptive quadrature, Laplace approximation, or Monte Carlo integration; for more complex nonlinear SDEs, particle-based likelihood methods may be needed. These approaches trade computational cost, approximation error, and implementation complexity. No approximation resolves structural non-identifiability: if the observation design cannot separate components, a more elaborate optimizer will not create the missing information.

Finally, individual random-effect estimates are conditional summaries, commonly posterior means or modes under the fitted population model. They are subject to shrinkage toward the population distribution, especially when a subject has few observations. Estimation of population parameters and prediction of individual effects are connected but distinct tasks.

Key takeaway: A mixed-effects SDE fit integrates each subject’s joint observation density over that subject’s shared random effect. Separate the random-effect, process, and measurement contributions in the model, enforce valid variance parameterizations, and verify that the numerical approximation and data design support the estimates.

Check Your Understanding

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

In the log-concentration model, which term contributes off-diagonal covariance between observations at different times: measurement error or Brownian variation?
Why is the subject likelihood an integral over one random effect rather than a product of integrals, one for each observation?
With the covariance decomposition in the worked example, what would the estimated $\tau^2$ be if the first variance were $0.15$, while the second variance and covariance remained $0.17$ and $0.11$?
What does a subject-specific posterior mean of $\eta_i$ represent, and why is it not itself a population parameter estimate?

Solutions. (1) Brownian variation contributes $\sigma^2\min(t_j,t_\ell)$; independent measurement error contributes only to the diagonal. (2) One latent effect is shared across the subject’s observations, so their joint conditional density is integrated once over its population distribution. (3) The decomposition gives $\sigma^2=0.17-0.15=0.02$, $\omega^2=0.11-0.02=0.09$, and $\tau^2=0.15-0.09-0.02=0.04$. (4) It summarizes that subject’s effect conditional on the fitted population model and their observations; it is an individual-level estimate, whereas population parameters describe the distribution across subjects.