Here Be Dragons: Bimodal posteriors arise from numerical integration error in longitudinal models

TL;DR

This paper reveals that numerical integration errors in longitudinal differential equation models can cause bimodal posterior distributions, leading to non-uniqueness in parameter estimation, and proposes solutions to mitigate this issue.

Contribution

It identifies a new source of non-uniqueness in Bayesian parameter estimation caused by numerical integration errors and offers methods to detect and address bimodal posteriors.

Findings

Bimodal posteriors can arise from numerical integration errors in longitudinal models.

Theoretical and empirical evidence demonstrates the existence of multi-modality.

Markov Chain Monte Carlo methods with proper solutions can avoid bimodality.

Abstract

Longitudinal models with dynamics governed by differential equations may require numerical integration alongside parameter estimation. We have identified a situation where the numerical integration introduces error in such a way that it becomes a novel source of non-uniqueness in estimation. We obtain two very different sets of parameters, one of which is a good estimate of the true values and the other a very poor one. The two estimates have forward numerical projections statistically indistinguishable from each other because of numerical error. In such cases, the posterior distribution for parameters is bimodal, with a dominant mode closer to the true parameter value, and a second cluster around the errant value. We demonstrate that multi-modality exists both theoretically and empirically for an affine first order differential equation, that a simulation workflow can test for evidence of the issue more generally, and that Markov Chain Monte Carlo sampling with a suitable solution can avoid bimodality. The issue of multi-modal posteriors arising from numerical error has consequences for Bayesian inverse methods that rely on numerical integration more broadly.