Simulation-Based Calibration and Model Comparison · Tutorial 92 of 100

Comparing ODE, SDE, and Mixed-Effects Candidates

Learn to design simulations and use within-subject contrasts and between-subject variation to distinguish ODE, SDE, and mixed-effects PK models.

Advanced 11 min read

What You'll Learn

Define comparable ODE, SDE, mixed-effects, and combined candidate models.
Use repeated observations to distinguish process variability from random-effect heterogeneity.
Derive moment-based diagnostics for a simple log-concentration model.
Design simulation-based model comparisons and interpret their selection rates.
Recognize boundary, observation-error, and model-misspecification pitfalls.

Why Candidate Models Need Different Generative Tests

A deterministic ODE, an SDE, and a mixed-effects model can all produce concentration profiles that vary across observed subjects. That resemblance is not enough to decide which model describes the source of variation. An SDE represents within-subject process variability; random effects represent between-subject heterogeneity; measurement error represents neither. A useful comparison simulates data under mechanisms that separate these sources, then checks whether fitting the candidate models recovers the mechanism that generated the data.

The aim is not to make one model win on every simulated dataset. It is to learn which patterns each model can reproduce, how often a comparison selects the wrong candidate, and whether the data-collection schedule can distinguish the candidates at all. The definitions below specify a deliberately simple setting in which the distinction can be derived exactly.

Definition: Let $X_i(t)$ be the log concentration for subject $i$, observed at times $t_j$. Consider the model $$ X_i(t)=m(t;\theta)+b_i+\sigma W_i(t), $$ with observations $Y_{ij}=X_i(t_j)+\epsilon_{ij}$. Here $m(t;\theta)$ is the common deterministic trajectory, $b_i\sim N(0,\tau^2)$ is a subject-level random intercept, $W_i(t)$ is standard Brownian motion independent across subjects and independent of $b_i$, and $\epsilon_{ij}\sim N(0,r^2)$ is independent measurement error. The four candidates are: ODE ($\tau^2=\sigma^2=0$); SDE only ($\tau^2=0,\sigma^2>0$); mixed-effects ODE ($\tau^2>0,\sigma^2=0$); and mixed-effects SDE ($\tau^2>0,\sigma^2>0$). All candidates use the same observation-error model unless the comparison explicitly tests that assumption.

This is a comparison framework, not a claim that every PK model has additive Brownian noise on log concentration. The log scale avoids negative concentrations after exponentiation, while the random-intercept and Brownian terms provide a transparent way to isolate the two sources of latent variability. The deterministic trajectory and observation model must be appropriate to the application; otherwise a comparison may select the least-wrong candidate rather than an adequate one.

What Repeated Measurements Reveal

For a short derivation, take $m(t;\theta)$ as known and subtract it from each observation. At two times $t_1<t_2$, write $Z_{i1}=Y_{i1}-m(t_1;\theta)$, $Z_{i2}=Y_{i2}-m(t_2;\theta)$, and $D_i=Z_{i2}-Z_{i1}$. The random intercept cancels in the difference. Brownian increments after $t_1$ remain, as does the difference of the two measurement errors.

Result: Under the model above, with independent errors of variance $r^2$, $$ \operatorname{Var}(D_i)=\sigma^2(t_2-t_1)+2r^2, $$ and $$ \operatorname{Var}(Z_{i1})=\tau^2+\sigma^2t_1+r^2. $$ Thus the within-subject change informs process variance after accounting for measurement error, while the level at the first time contains both random-intercept and accumulated process variance.

Derivation. Since $W_i(t_2)-W_i(t_1)$ is independent of $W_i(t_1)$ and has variance $t_2-t_1$, $$ D_i=\sigma\{W_i(t_2)-W_i(t_1)\}+\epsilon_{i2}-\epsilon_{i1}. $$ The independent terms have variances $\sigma^2(t_2-t_1)$, $r^2$, and $r^2$, giving the first expression. Also, $$ Z_{i1}=b_i+\sigma W_i(t_1)+\epsilon_{i1}. $$ The three terms are independent and have variances $\tau^2$, $\sigma^2t_1$, and $r^2$, giving the second expression. These equations show why repeated observations at separated times matter: a between-subject intercept disappears from a within-subject difference, but process noise accumulates over the interval.

The derivation also gives a simple identification check. If $r^2$ is known and $t_1>0$, the two population variances determine the two variance components: $$ \sigma^2=\frac{\operatorname{Var}(D_i)-2r^2}{t_2-t_1}, \qquad \tau^2=\operatorname{Var}(Z_{i1})-\sigma^2t_1-r^2. $$ A negative estimate of either variance component is not evidence for a physically negative variance. It signals sampling variation, an incompatible model, or an inadequate error assumption; a constrained fit places that component at zero.

Worked Example: Separating Process Noise From a Random Intercept

Suppose log concentrations are measured at $t_1=1$ and $t_2=2$, with known measurement-error variance $r^2=0.04$. After subtracting the fitted deterministic trajectory, a simulated dataset has sample variance $S_D^2=0.44$ for the paired differences and $S_{Z_1}^2=0.68$ for first-time residuals. Use the moment equations to estimate the two components.

The interval is $t_2-t_1=1$, so the process-variance estimate is $$ \hat{\sigma}^2=\frac{S_D^2-2r^2}{t_2-t_1} =\frac{0.44-2(0.04)}{1} =0.36. $$ The two measurement errors contribute $0.08$ to the variance of a difference, leaving $0.36$ for process variation across a one-time-unit interval.

Next, $$ \hat{\tau}^2=S_{Z_1}^2-\hat{\sigma}^2t_1-r^2 =0.68-(0.36)(1)-0.04 =0.28. $$ Both estimates are nonnegative. Under the fitted combined model, the predicted first-time variance is $0.28+0.36+0.04=0.68$, matching the supplied summary. The predicted difference variance is $0.36+0.08=0.44$, also matching it. An SDE-only candidate would omit the $0.28$ random-intercept contribution; a mixed-effects ODE would predict only $0.08$ for the difference variance, the measurement-error contribution.

Moment Estimates and Their Assumptions

The sample-variance calculation is a diagnostic and a transparent starting point, not a replacement for fitting the full model. Its validity depends on independent subjects, the specified covariance structure, a correctly removed mean trajectory, and known $r^2$. If $m(t;\theta)$ is estimated from the same data, uncertainty in that estimate can affect the residual summaries. If observation error is unknown, process noise and measurement error can be hard to distinguish, especially when observations are closely spaced.

Result: Suppose subjects are independent and identically distributed under the two-time model, $r^2$ is known, and $t_2>t_1$. Define the unbiased sample variances $S_D^2$ and $S_{Z_1}^2$ using denominator $n-1$. Then $$ \tilde{\sigma}^2=\frac{S_D^2-2r^2}{t_2-t_1}, \qquad \tilde{\tau}^2=S_{Z_1}^2-\tilde{\sigma}^2t_1-r^2 $$ are unbiased for $\sigma^2$ and $\tau^2$, respectively. They need not be nonnegative in a finite sample.

Proof. By the variance expressions above, the sample variance is unbiased for the population variance, so $$ E[S_D^2]=\sigma^2(t_2-t_1)+2r^2. $$ Linearity of expectation gives $$ E[\tilde{\sigma}^2] =\frac{E[S_D^2]-2r^2}{t_2-t_1} =\sigma^2. $$ Likewise, $E[S_{Z_1}^2]=\tau^2+\sigma^2t_1+r^2$. Therefore $$ E[\tilde{\tau}^2] =E[S_{Z_1}^2]-t_1E[\tilde{\sigma}^2]-r^2 =\tau^2+\sigma^2t_1+r^2-\sigma^2t_1-r^2 =\tau^2. $$ Unbiasedness does not impose the nonnegative parameter constraints of a variance-component model. A constrained likelihood or other constrained estimator may therefore differ from these moment estimates. $\square$

Worked Example: A Pattern Consistent With Between-Subject Heterogeneity

Keep $t_1=1$, $t_2=2$, and $r^2=0.04$. Now suppose a different simulated dataset gives $S_D^2=0.08$ and $S_{Z_1}^2=0.68$. Calculate the components and interpret the pattern.

The difference variance is entirely accounted for by the two observation errors: $$ \tilde{\sigma}^2=\frac{0.08-2(0.04)}{1}=0. $$ Then $$ \tilde{\tau}^2=0.68-(0)(1)-0.04=0.64. $$ The summary is consistent with substantial between-subject heterogeneity and no detectable process variance: subjects differ in their levels, but the paired changes show no extra variability beyond measurement error. A mixed-effects ODE can represent this pattern. An ODE with no random effects predicts no between-subject latent variation, while an SDE-only model attributes variation to within-subject changes and does not explain the large first-time variance through a random intercept. These summaries are not a formal selection test; sampling uncertainty and the full time-course likelihood still matter.

Designing a Simulation-Based Comparison

A fair comparison begins with the same observation schedule, dosing information, residual-error model, and estimation protocol for every candidate. Generate separate simulation scenarios under the ODE, SDE-only, mixed-effects ODE, and combined model. Include settings with only one source of latent variation, settings with both, and values near zero where separation is difficult. For each simulated dataset, fit every candidate using its appropriate likelihood, integrating over random effects and latent paths as required. Compare models using a prespecified criterion, such as a common marginal likelihood-based information criterion or held-out predictive log score.

Specify generating mechanisms.
Set the structural trajectory, observation times, observation error, sample size, and variance components. Vary process and random-effect variances separately as well as jointly.

Simulate at the subject level.
For each subject, draw one random effect, then a process path at the scheduled times, then observation errors. Keeping these draws distinct prevents accidental substitution of measurement noise for process noise.

Refit all candidates consistently.
Use the same data and mean-trajectory treatment for each fit. Record convergence, boundary estimates, selected model, and predictive or likelihood-based comparison scores.

Summarize recovery and failure.
For each generating scenario, tabulate selection frequencies and parameter-estimation behavior. A useful candidate should not only fit its own scenario but should also avoid frequent false selection under competing mechanisms.

Worked Example: Interpreting Selection Frequencies

Consider a planned comparison with 1,000 simulated datasets generated under each of four scenarios. Suppose the fitted-model selection frequencies in the process-only scenario are: ODE 0.06, SDE-only 0.71, mixed-effects ODE 0.07, and combined model 0.16. These are illustrative frequencies for interpreting a simulation output, not a universal benchmark.

The process-only generating candidate is selected most often, so the design contains information about process variability. However, the SDE-only model is selected in only 71% of replicates; a 16% selection rate for the combined model suggests that finite-sample data sometimes support an unnecessary random effect. Under this scenario, the estimated Monte Carlo standard error for the 0.16 frequency is $$ \sqrt{\frac{0.16(1-0.16)}{1000}}\approx 0.0116. $$ Thus simulation uncertainty alone is about 1.2 percentage points. The observed 16% is not plausibly explained by that small Monte Carlo error if the true frequency were near zero; it is a meaningful property of the fitting-and-selection procedure that deserves investigation.

Repeat the calculation for data generated under each candidate. In particular, a combined model selected frequently when $\tau^2=0$ or $\sigma^2=0$ may be overfitting a variance component. Also inspect parameter estimates and predictive checks: the correct model can be selected while its variance estimate remains biased or highly variable.

Comparison, Diagnostics, and Limits

Candidate	Latent source of variation	Characteristic two-time pattern	Typical diagnostic concern
ODE	None; only the shared deterministic trajectory	Variation explained by observation error	Repeated structured deviations exceed the error model
SDE only	Within-subject process noise	Excess variance in changes, increasing with interval length	Persistent subject-level offsets remain unexplained
Mixed-effects ODE	Between-subject random effects	Subject levels vary; changes have no extra process variance	Within-subject changes vary beyond observation error
Mixed-effects SDE	Both sources	Heterogeneous levels and excess within-subject change variance	Extra flexibility may be weakly identified or unnecessary

A likelihood-ratio comparison of a model with a variance component against one where that component is zero requires special care: zero lies on the boundary of the allowable parameter space. The usual chi-square reference distribution for a regular, interior null parameter need not apply. A parametric bootstrap that simulates under the restricted model and repeats the fitting and comparison procedure can calibrate the test in the setting actually used. Information criteria and predictive scores also require consistent likelihood definitions; comparing a conditional likelihood for one model with a marginal likelihood for another is not an apples-to-apples comparison.

Diagnostics should be tied to the simulated mechanism. Check whether within-subject innovations retain serial structure, whether subject-level residual summaries show unexplained dispersion, and whether fitted variance components pile up at zero or change markedly with small modeling choices. Check convergence and path or integration accuracy as well: a numerical fitting failure is not evidence that a candidate is scientifically wrong. The earlier material on simulation calibration provides the general logic for assessing whether a procedure behaves as intended.

Common Pitfalls

Calling all unexplained spread “process noise.” Persistent subject-to-subject differences can inflate pooled residual variation. Use repeated measurements and within-subject contrasts to test whether variation occurs over time within subjects.
Treating measurement error as process variability. Independent measurement errors contribute $2r^2$ to a paired difference. Ignoring that contribution overstates the process variance.
Assuming a single schedule distinguishes every candidate. Closely spaced observations may contain little information about gradual process variation, while sparse observations may conceal its temporal pattern. Simulate the actual schedule.
Choosing the combined model because it is most flexible. A variance component estimated near zero can add complexity without improving meaningful prediction. Evaluate selection under scenarios where that component is truly absent.
Interpreting one simulation replicate as a verdict. Model selection is variable across datasets. Use repeated simulations and report selection frequencies and their Monte Carlo uncertainty.

The main limitation of this comparison is that the simple random-intercept and Brownian model encodes a particular covariance pattern. Random slopes, time-varying covariates, misspecified structural kinetics, correlated measurement errors, or non-Brownian process noise can imitate or obscure the patterns described here. The simulation scenarios should therefore reflect plausible alternatives, not only the four idealized candidates. The purpose is to assess discrimination under stated assumptions, not to prove that one stochastic mechanism is uniquely true.

Key takeaway: Compare ODE, SDE, and mixed-effects candidates by simulating the mechanisms separately and examining both within-subject changes and between-subject levels. A defensible conclusion requires a shared observation model, calibrated selection behavior, and diagnostics that match the proposed source of variability.

Check Your Understanding

Answer each question, then compare with the solutions.

With $t_1=1$, $t_2=3$, $r^2=0.02$, and $\operatorname{Var}(D)=0.28$, what is the implied process variance?
Why can a random intercept not, by itself, increase the variance of a within-subject difference?
In a simulation under a mixed-effects ODE, what pattern would indicate that an SDE component is being spuriously selected?
Why is the usual chi-square likelihood-ratio reference potentially unsuitable for testing whether a variance component is zero?

Solutions. (1) The interval is $2$, and the difference contains measurement-error variance $2r^2=0.04$. Therefore $\sigma^2=(0.28-0.04)/2=0.12$. (2) The same subject intercept appears in both observations and cancels: $(b_i+\cdots)-(b_i+\cdots)$ contains no $b_i$. (3) Frequent selection of the combined model, or positive estimated process variance beyond the observation-error contribution, indicates possible spurious selection; assess its frequency across replicates rather than from one dataset. (4) The null value zero is on the boundary of the nonnegative variance-parameter space, violating the regular interior-parameter conditions for the usual chi-square result; simulation-based calibration can assess the reference behavior for the actual procedure.

Back to Tutorials