Bayesian inference for ordinary differential equations models with heteroscedastic measurement error

1. Problem Statement

Ordinary Differential Equation (ODE) models are fundamental for describing dynamical systems in biology, ecology, and epidemiology. A critical step in utilizing these models is parameter estimation from observational data. Standard Bayesian inference for ODEs typically assumes that measurement errors are homoscedastic (constant variance) and independent, following a Gaussian distribution.

However, real-world data often exhibit heteroscedasticity, where the variance of measurement error changes over time or depends on the system state. This can arise from:

• Environmental fluctuations.
• Changes in measurement instrument resolution.
• Intrinsic system variability (e.g., higher variance during epidemic peaks).
• Model misspecification.

The Core Issue: Assuming homoscedasticity when heteroscedasticity exists leads to:

• Biased parameter estimates.
• Inaccurate uncertainty quantification (e.g., overconfident posterior distributions that are too narrow).
• Poor predictive performance, particularly in regions of high or low measurement variability.

While parametric heteroscedastic models exist, they require specifying a functional form for the variance, which risks misspecification and increases computational cost due to additional parameters to estimate within the Markov Chain Monte Carlo (MCMC) loop.

2. Methodology: A Two-Step Semi-Parametric Framework

The authors propose a novel two-step semi-parametric framework that decouples the estimation of the error structure from the inference of the ODE parameters. This approach avoids the computational burden of joint estimation while maintaining flexibility.

Step 1: Estimation of the Measurement Error Process (HetGP)

• Goal: Learn the time-dependent error structure directly from the data without assuming a specific parametric form.
• Technique: The authors employ Heteroscedastic Gaussian Processes (HetGP).
  • The latent mean function $f(t)$ and the log-variance function $\log \sigma^2(t)$ are modeled as two coupled Gaussian processes.
  • The log-transformation ensures the variance remains positive.
• Algorithm: An alternating optimization scheme is used to maximize the coupled log-marginal likelihoods:
  1. Estimate the mean GP hyperparameters given a current variance estimate.
  2. Compute squared residuals between observations and the posterior mean.
  3. Fit a second GP to the log-squared residuals to update the variance hyperparameters.
  4. Iterate until convergence.
• Output: A time-dependent variance estimate $\hat{\sigma}^2(t_i)$ for each observation time. The framework includes an automatic model selection feature (checkHom) to compare the marginal likelihood of the heteroscedastic model against a homoscedastic one.

Step 2: Bayesian Inference for ODE Parameters

• Goal: Estimate the ODE parameters $\theta$ using the learned error structure.
• Technique: The estimated variance $\hat{\sigma}^2(t_i)$ from Step 1 is treated as known and incorporated into the likelihood function.
• Likelihood: A heteroscedastic Gaussian likelihood is constructed:
  $$\log p(\tilde{y} | \theta, \hat{\sigma}^2(t)) = -\frac{1}{2} \sum_{i=1}^N \left[ \log(2\pi\hat{\sigma}^2(t_i)) + \frac{(\tilde{y}(t_i) - g(y(t_i; \theta)))^2}{\hat{\sigma}^2(t_i)} \right]$$
• Inference: Standard Metropolis-Hastings (MH) MCMC is used to sample from the posterior distribution $p(\theta | \tilde{y})$. Because the error structure is pre-estimated, the MCMC sampler does not need to jointly optimize the error parameters, significantly improving computational efficiency and stability.

3. Key Contributions

1. Decoupled Inference Strategy: The primary contribution is the separation of stochastic error modeling (via non-parametric HetGP) from mechanistic parameter inference. This avoids the "curse of dimensionality" and identifiability issues often found in joint estimation.
2. Semi-Parametric Flexibility: The method does not require the user to specify a parametric form for the variance (e.g., linear or exponential dependence on time). The error structure is learned directly from the data.
3. Automated Model Selection: The integration of the checkHom argument allows the framework to automatically determine if a heteroscedastic model is statistically preferred over a homoscedastic one based on marginal likelihood.
4. Robust Uncertainty Quantification: The framework provides predictive intervals that dynamically adapt to the data's variability, avoiding the overconfidence associated with constant-variance assumptions.

4. Results

The authors validated the framework through three case studies:

• Case 1: Logistic Growth Simulation (Synthetic Data)

  • Setup: Simulated logistic growth with known time-varying error.
  • Metric: Maximum Mean Discrepancy (MMD) between posterior samples and the true posterior.
  • Finding: As the number of observations increased, the homoscedastic model's posterior diverged significantly from the truth (biased). The proposed HetGP approach consistently produced posteriors closer to the true distribution, with MMD decreasing as sample size grew.

• Case 2: Coral Reef Recovery (Real Ecological Data)

  • Setup: Modeling coral cover recovery using a Richards' growth model.
  • Finding: The HetGP identified systematic changes in measurement variability over time. While marginal parameter distributions ($\alpha, \gamma, K$) appeared similar between models, the joint posterior distributions showed different correlation structures. Crucially, the predictive intervals from the proposed method expanded and contracted to match the data's variability, whereas the homoscedastic model produced uniform, often misleading, intervals.

• Case 3: 1984 Measles Epidemic (Real Epidemiological Data)

  • Setup: SIR model for measles incidence in England and Wales.
  • Finding: The error variance was significantly higher during epidemic peaks. The proposed method yielded more concentrated posterior distributions for the transmission rate ($\beta$) and initial susceptible population ($S_0$).
  • Impact: The basic reproduction number ($R_0$) estimated by the proposed method was more concentrated and lower than the homoscedastic estimate, aligning better with the observed epidemic timing and scale. The homoscedastic model produced overly broad and shifted $R_0$ distributions.

5. Significance and Conclusion

This work demonstrates that ignoring heteroscedasticity is not a minor technicality but a source of significant bias in ODE modeling. The proposed two-step framework offers a practical, automated, and computationally efficient solution for handling time-varying measurement errors.

• Scientific Impact: It enables more reliable parameter estimation and uncertainty quantification, which is critical for decision-making in fields like public health (e.g., predicting epidemic thresholds) and ecology (e.g., forecasting species recovery).
• Future Directions: The authors note that while the current implementation assumes Gaussian errors, the framework is general enough to accommodate non-Gaussian error structures (e.g., via Generalized Additive Models) in future work to handle skewed or heavy-tailed data.

In summary, the paper provides a robust workflow for Bayesian ODE inference that bridges the gap between flexible statistical learning and mechanistic modeling, ensuring that uncertainty is quantified accurately in the presence of complex, real-world noise.