Why There Are Two Integrals
In a mixed-effects SDE model, each subject has a random effect that is not observed, and a latent state that evolves stochastically between observations. Even if the population parameters are fixed, neither the subject effect nor the state path is known. The likelihood for that subject must account for both sources of uncertainty. Keeping the two integrations distinct is essential: a method that handles the latent path well does not automatically handle the random-effects integral, and vice versa.
Let \(Y_i=(Y_{i1},\ldots,Y_{im_i})\) denote subject \(i\)'s observations at times \(t_{i1},\ldots,t_{im_i}\). Let \(b_i\) be that subject's random-effect vector, with density \(g_\theta(b_i)\). Conditional on \(b_i\), suppose the state \(X_i(t)\) follows the specified SDE model, whose probability law on paths is denoted \(P_\theta^{b_i}\). The parameter vector \(\theta\) includes population, observation, and SDE parameters. Assume observations are conditionally independent given the state at their observation times and \(b_i\), with density \(q_\theta(y\mid x,b_i)\).
The notation \(P_\theta^b(dX)\) is deliberate. An SDE path is generally not represented by an ordinary density with respect to a useful infinite-dimensional Lebesgue measure. The inner integral is an expectation under a probability law on paths. In a discretized implementation, one may instead integrate over a finite vector of states, but that is an approximation or a model defined on the discretized grid.
What Can Be Integrated, and in What Order?
The integrand in the likelihood is nonnegative, and Tonelli's theorem applies to the iterated integral defined by the joint law of \(b\) and \(X\). To write that integral in the reverse order, use a disintegration of the joint law: the inner measure is the conditional law of \(b\) given \(X\), not generally its marginal density \(g_\theta(b)\,db\). For a proper observation model, the marginal density is finite almost everywhere under the model. This fact permits a rearrangement; it does not make either integral easy to calculate.
Proof. By the definition of the joint law \(P_\theta\), the first integral is \(\int H(b,X)\,P_\theta(db,dX)\). Disintegrating this joint law with respect to its path marginal gives \(\int P_{\theta,X}(dX)\int H(b,X)\,P_\theta(db\mid X)\). The conditional law \(P_\theta(db\mid X)\) accounts for the dependence between the effect and path. \(\square\)
For computation, it is often helpful to separate the path integral into operations at observation times. Write \(K_{\theta,j}^b(dx_j\mid x_{j-1})\) for the SDE transition probability from the preceding observation time to \(t_{ij}\), conditional on \(b\). It includes the evolution over the interval between observations. If \(x_0\) has conditional distribution \(\mu_\theta^b\), the conditional likelihood is
This identity follows by applying the Markov property to the state at successive observation times and multiplying each transition measure by the observation density at its endpoint. Integrating over all the \(x_j\)'s sums over every possible sequence of observation-time states. The unobserved within-interval portions of the path have already been integrated into the transition probabilities. Known doses or inputs must be included in the appropriate transition kernel; they do not remove either marginalization.
Worked Example: Exact Gaussian Marginalization
Let \(B\sim N(0,1)\) be a subject effect. Conditional on \(B=b\), let the latent endpoint \(X\) have distribution \(N(b,2)\), and suppose the observed value is \(Y=X\). Find the marginal density of \(Y\) at \(y=1\).
Inner integration. Given \(b\), the observation density is $$ p(y=1\mid b)=\frac{1}{\sqrt{4\pi}}\exp\left[-\frac{(1-b)^2}{4}\right]. $$
Outer integration. Since \(Y=B+\varepsilon\), where \(\varepsilon\sim N(0,2)\) is independent of \(B\), the sum is \(Y\sim N(0,1+2)=N(0,3)\). Thus the nested integral is exactly $$ p(Y=1) = \int \frac{1}{\sqrt{2\pi}}e^{-b^2/2} \frac{1}{\sqrt{4\pi}}e^{-(1-b)^2/4}\,db = \frac{1}{\sqrt{6\pi}}e^{-1/6} \approx 0.1950. $$
The variance is \(3\), not \(2\): the random effect contributes variance \(1\), and the conditional latent endpoint contributes variance \(2\). The calculation illustrates both levels of the likelihood, even though Gaussian closure lets us perform them analytically.
How Estimation Strategies Approximate the Nesting
When exact integration is unavailable, an estimation strategy must say what it does with each integral. In the table, “inner” means the path integral conditional on a specified \(b\), and “outer” means averaging over \(g_\theta(b)\). An algorithm may approximate both, or use a stochastic calculation at one level while treating the other differently.
| Strategy | Path integral, conditional on \(b\) | Random-effects integral | What to check |
|---|---|---|---|
| Numerical transition integration | Integrates or approximates successive state transitions | Quadrature or Monte Carlo over \(b\) | Grid and quadrature accuracy at both levels |
| Laplace approximation | May use a joint approximation over selected latent variables | Approximates an integral near a conditional or joint mode | Whether the integrand is sufficiently concentrated and locally quadratic |
| Particle filter inside outer integration | Uses particles to estimate the conditional likelihood for fixed \(b\) | Quadrature or sampling over \(b\) | Particle variability as well as outer integration error |
| Particle MCMC or related Bayesian methods | Uses particles or latent-state simulation in posterior computation | Samples effects and population parameters jointly or conditionally | Monte Carlo mixing and whether the target includes both latent layers |
| Conditional-mode approximations | Often simplifies or approximates state uncertainty | Approximates effect integration around a mode | Bias from both the mode approximation and any state simplification |
A particle filter run at fixed \(b\) targets the inner conditional state calculation. It does not, by itself, integrate over the population distribution of \(b\). A Bayesian sampler that draws both \(b\) and \(X\) can target their joint posterior, but posterior simulation is not the same numerical task as evaluating the marginal likelihood. Likewise, a Laplace approximation is not an exact integration: its accuracy depends on how well a local quadratic expansion represents the relevant integrand.
Worked Example: Three-Point Quadrature for an Outer Integral
Suppose that after the latent-path calculation, the conditional likelihood is proportional to \(h(b)=\exp(-b^2/2)\), and \(B\sim N(0,1)\). The required outer expectation is \(E[h(B)]\). A three-node Gauss-Hermite rule for a standard normal expectation uses nodes \(0,\sqrt{3},-\sqrt{3}\), with weights \(2/3,1/6,1/6\).
Applying the rule gives $$ \widehat E[h(B)] = \frac{2}{3}h(0)+\frac{1}{6}h(\sqrt{3})+\frac{1}{6}h(-\sqrt{3}) = \frac{2}{3}+\frac{1}{3}e^{-3/2} \approx 0.741043. $$
The exact value follows by combining the standard normal density with \(h\): $$ E[h(B)] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-b^2}\,db = \frac{1}{\sqrt{2}} \approx 0.707107. $$
The absolute error for this quadrature rule is therefore $$ \frac{2}{3}+\frac{e^{-3/2}}{3}-\frac{1}{\sqrt{2}} \approx 0.033937. $$ This is error in the outer integral in this example; it says nothing about the accuracy of a separate approximation used to calculate the inner path contribution.
Worked Example: Inner Sum Followed by Outer Average
A small finite-state model makes the order of operations visible. Let the random effect take values \(b=+1\) and \(b=-1\), each with probability \(1/2\). A latent state \(X\) is binary, with \(P(X=1\mid b=+1)=0.8\) and \(P(X=1\mid b=-1)=0.2\). The observation is binary with \(P(Y=1\mid X=1)=0.9\) and \(P(Y=1\mid X=0)=0.1\). Calculate the likelihood of observing \(Y=1\).
For \(b=+1\), integrate over the two possible latent states: $$ P(Y=1\mid b=+1)=0.8(0.9)+0.2(0.1)=0.72+0.02=0.74. $$ For \(b=-1\): $$ P(Y=1\mid b=-1)=0.2(0.9)+0.8(0.1)=0.18+0.08=0.26. $$
Now average the conditional likelihood over the random effects: $$ P(Y=1)=\frac{1}{2}(0.74)+\frac{1}{2}(0.26)=0.50. $$ The calculation first sums over the latent state conditional on each effect, then averages over effects. The joint alternatives give the same result: the probabilities for \((b,X)=(+1,1),(+1,0),(-1,1),(-1,0)\) are \(0.4,0.1,0.1,0.4\), whose contributions to \(P(Y=1)\) are \(0.36,0.01,0.09,0.04\), totaling \(0.50\).
Diagnostics, Limitations, and a Common Pitfall
A useful diagnostic is to vary the numerical resolution at each level separately. Increase outer quadrature nodes while holding the inner calculation fixed; then increase particles or refine the state grid while holding the outer rule fixed. If the likelihood or parameter estimates move substantially in either experiment, the corresponding integration is not yet stable. For particle-based inner calculations, repeat runs at fixed \(b\) to assess Monte Carlo variability. For Laplace approximations, inspect whether the conditional integrand is sharply concentrated and approximately symmetric around its mode; skewness, multiple modes, or a broad tail can make a local approximation unreliable.
A common mistake is to treat a conditional likelihood at a typical effect, such as \(b=0\), as the marginal likelihood. In general, $$ L_i(\theta)\ne L_i(\theta\mid b=0), $$ because an average of likelihoods over the random-effects distribution is not generally equal to the likelihood at the average effect. The correction is to perform the outer integral, exactly or by a stated approximation. A second pitfall is to claim that a particle filter has “integrated the random effects” merely because it represents uncertainty in the state: unless the effect is included in the particle state or separately integrated, the filter calculation is conditional on the chosen effect.
No integration method can repair a misspecified random-effects distribution, observation model, or SDE. Numerical stability is evidence about computation, not proof that the model is appropriate. The practical comparison is therefore between explicit approximations: quadrature can be effective for low-dimensional effects; particle methods are useful for nonlinear state dynamics but introduce Monte Carlo error; and Laplace methods can be efficient when concentration near a well-behaved mode is credible. High-dimensional random effects or difficult latent dynamics may make any one method inadequate without further approximation or computation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- If a particle filter is run at one fixed value \(b=b^\ast\), which part of the nested likelihood does it address?
- In the Gaussian example, what is the marginal variance of \(Y\), and why is it not the conditional variance?
- For the three-point quadrature example, identify the exact outer expectation and the absolute quadrature error to six decimal places.
- In the binary-state example, what is \(P(Y=1\mid b=+1)\), and what is the fully marginalized \(P(Y=1)\)?
Solutions. (1) It estimates or approximates the inner path integral conditional on \(b^\ast\); an additional calculation is needed to integrate over \(g_\theta(b)\). (2) The marginal variance is \(1+2=3\), because \(Y=B+\varepsilon\), with independent variances \(1\) and \(2\); the conditional variance given \(B\) is only \(2\). (3) The exact value is \(1/\sqrt{2}\), and the absolute error is \(2/3+e^{-3/2}/3-1/\sqrt{2}\approx 0.033937\). (4) The conditional probability is \(0.74\), and averaging it with the \(b=-1\) value \(0.26\), using equal weights, gives \(P(Y=1)=0.50\).