Using flexible noise models to avoid noise model misspecification in inference of differential equation time series models

The dataset consists of time points $\{t_{i}\}_{i=1}^{N}$ and corresponding noisy data $\{y_{i}\}_{i=1}^{N}$. A typical modelling assumption is to treat the noise on each data point as IID Gaussian with a variance parameter $\sigma^{2}$, so that,

Our first step is to generalize eq. (3) so that the variance of noise terms can vary (i.e. allow the noise to be non-identically distributed), and each noise realisation can be correlated with its neighbours (i.e. be non-independent). A multivariate Gaussian can handle both of these generalisations, where we model a random vector, $\mathbf{y}=(y_{1},\dots,y_{N})^{\top}$, as having a mean, $\mathbf{f}(\theta)=(f(t_{1};\theta),\dots,f(t_{N};\theta))^{\top}$,

For appropriate values of the covariance matrix $\Sigma$, this distributional assumption encompasses a wide variety of noise forms which may include correlated and heteroscedastic noise terms. For example, eq. (4) could describe heteroscedastic noise which scales with the magnitude of the trajectory with $\Sigma=\text{diag}(\mathbf{f}(\theta)\sigma^{2})$. For autocorrelated noise terms, $\Sigma$ would be dense.

Multiple methods have been proposed for inference of covariance matrices (Hoffbeck and Landgrebe, 1996; Lam and Fan, 2009; Diggle and Verbyla, 1998; Bickel and Levina, 2008; Cai and Liu, 2011; Schäfer and Strimmer, 2005). A standard Bayesian approach places a prior on $\Sigma$ and infers it along with ODE model parameters, $\theta$. Typical choices for priors include the conjugate inverse-Wishart (Gelman et al., 2013; Huang and Wand, 2013), or a prior based around the LKJ correlation matrix (Lewandowski et al., 2009; Stan Development Team, 2016). However, these methods are not designed to handle the covariance of a single time series. For a single time series obeying eq. (4), there is just one multivariate data point (that is, the vector $\mathbf{y}$) available to inform the matrix $\Sigma$. With such limited data, standard methods for estimating covariance matrices have too much freedom, resulting in dense matrices that overfit the data.

Our strategy is to impose a positive definite covariance function $C:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$, which generates a covariance matrix according to the rule,

For example, heteroscedastic errors, where $\Sigma=\text{diag}(\mathbf{f}(\theta)\sigma^{2})$, could be represented by the following covariance function:

where $\delta_{ij}=1,\;\text{if }i=j;\;0,\text{otherwise}$. In this paper, we consider positive definite kernels which are flexible enough to capture a wide variety of noise forms, with parameters that can, nonetheless, be learned from a single time series.

In this section, we introduce the kernels used throughout this paper. Notwithstanding the important differences discussed in §2.4, much of the work on kernel functions for Gaussian processes is applicable to ODE noise models as well, and the three kernels we discuss have seen extensive use in Gaussian process regression. One of the most widely used positive definite kernels is the Gaussian kernel (also called the Radial Basis Function, or RBF) (Fasshauer, 2011),

We also consider the Laplacian kernel (Feragen et al., 2015) for specifying time series autocovariances, since it more faithfully reproduces the types of persistence emergent from basic univariate time series models,

The kernels in eqs. (7) & (8) are each characterised by two parameters which control the size and autocovariance in the errors. A more general class of kernels is the Matérn (Williams and Rasmussen, 2006),

where $K_{\nu}$ is the modified Bessel function of the second kind. For $\nu=1/2$, the Matérn kernel simplifies to the Laplacian kernel.

Consider a function $g:\mathcal{X}\rightarrow\mathbb{R}$ obeying a GP with mean function $m$ and kernel $C$, i.e. $g\sim\mathcal{GP}(m,C)$ (see, for example, Rasmussen (2003)). For every finite set of inputs $\{t_{i}\}_{i=1}^{N},\,t_{i}\in\mathcal{X}$, the vector of function values $\mathbf{g}=(g(t_{1}),\dots,g(t_{N}))^{\top}$ has a multivariate Gaussian distribution,

where $\mathbf{m}=(m(t_{1}),\dots,m(t_{N}))^{\top}$ and $\Sigma$ is generated as in eq. (5). This distribution, identical with eq. (4) for $m(\cdot)=f(\,\cdot\,;\theta)$, illustrates an apparent resemblance between the multivariate normal likelihood for time series noise and the GP. Our proposed noise model, however, differs from a GP in several key aspects:

1. 1.

   In GP regression, eq. (10) determines a prior over functions, and the posterior over functions is inferred. Our proposed noise model uses the multivariate normal specification as a likelihood for finite observed data, and posterior inference applies only to the parameters of $f$, not the functional form of the noise-free relationship between $y$ and $t$ which is assumed fixed and fully determined by $\theta$.

2. 2.

   To handle noisy data, Gaussian process regression typically adds an extra noise term—often IID Gaussian. No such terms are used in our multivariate normal noise process.

That is, in full, the likelihood for our multivariate normal model is given by eq. (4), with covariance matrix given by eq. (5). An example of the utility of the multivariate normal noise process is shown in Figure S1. In this example, we show that the Laplacian kernel can faithfully capture autoregressive order 1 (AR(1)) noise in an ODE time series model, enabling accurate posterior inference for the ODE model parameters.

Non-stationary covariance functions have been used for spatial modelling and Gaussian process regression. Unlike stationary kernels which depend only on the distance between the two inputs, in the non-stationary case, the kernel shape itself must depend on the input location. This is expressed using the notation $k_{s}(u)$ for a kernel centred at location $s$ and evaluated at location $u$. For example, for the Laplacian kernel with one dimensional input $t$, we would take

with the kernel parameters $\sigma$ and $L$ being functions of the kernel centre location $t$.

If eq. (11) is used to construct a covariance matrix, there are no guarantees that it will be positive definite. Instead, non-stationary modelling has relied on the following general formula for a non-stationary positive definite covariance function:

for inputs $x_{i},x_{j},u\in\mathbb{R}^{N}$ (Higdon et al., 1999; Paciorek, 2003). The non-stationary version of the Gaussian RBF covariance function can be derived from this formula, which has been used in non-stationary Gaussian process regression (Gibbs, 1998; Paciorek and Schervish, 2004). To learn time-varying kernel parameters, a Gaussian process prior can be placed on each (Paciorek and Schervish, 2004; Heinonen et al., 2016).

In this section, we present a non-stationary version of the Laplacian kernel. The techniques presented here are equally applicable to any appropriate positive definite kernel, however.

The one-dimensional non-stationary Laplacian covariance function is:

Eq. (13) may be derived as a special case of the non-stationary Matérn kernel (Paciorek and Schervish, 2004); it also follows directly from the one-dimensional case of eq. (12) using reparametrised versions of the respective stationary kernels, with the reparametrisations chosen to ensure that the final non-stationary covariance functions have a sensible form (cf. eq. (3.69) in Gibbs (1998)). The logarithms of $L$ and $\sigma$ each vary over time governed by Gaussian process priors:

where $\mu_{l}$, $K_{l}$, $\mu_{\sigma}$, and $K_{\sigma}$ are the GP hyperparameters.

Having specified a non-stationary covariance function such as eq. (13), the next task is to infer the posterior distribution of model and covariance parameters. However, analytic expressions for the posterior mean and variance of the Gaussian processes $L(t)$ and $\sigma(t)$ are not available. Instead, MCMC sampling or maximum a posteriori (MAP) estimation can be used to infer the values of $L(t)$ and $\sigma(t)$ at each time point (Heinonen et al., 2016). MCMC sampling yields a set of samples distributed according to the posterior distribution, while MAP estimation uses optimisation to find the parameters values at which the posterior distribution is maximised. For both MCMC and MAP estimation, we recommend the use of gradient-based methods (e.g., Hamiltonian MCMC and gradient-based optimisers) for improved convergence rates in the high dimensional parameter space (Neal, 2011). When analytic gradients are not available, automatic differentiation can be used. Indeed, all our GP examples presented here involve an interpolation scheme discussed in §5, rendering analytic derivatives intractable, and we resort to using automatic differentiation.

MCMC sampling is advantageous because it can, in theory, fully characterise posterior uncertainty in the covariance parameter fits. However, it is computationally costly and may be not be necessary for concentrated unimodel posteriors. In these cases, MAP estimation may be preferred due to its lower computational burden. While MAP estimates are often sufficient for the kernel parameter fits, we typically still require information about the posterior uncertainty in the ODE model parameters. Thus, we recommend the following procedure for long ODE time series problems with non-stationary covariance functions, which is specified in Algorithm 1. First, the joint MAP estimate of model parameters and covariance parameters is obtained using a gradient-based optimiser. Then, MCMC sampling is used to obtain the posterior distribution of model parameters conditional on the previously obtained MAP estimate of covariance parameters. For both optimisation and MCMC sampling, random or uniform initialisation of the covariance parameters will work for easier problems but will delay convergence on longer time series. In long time series problems with intelligible noise patterns, we recommend a data-driven initialisation of the covariance parameters in order to accelerate convergence of MCMC or optimisation. To initialise $L$ and $\sigma$ in the non-stationary Laplacian covariance function, we use the procedure given by Algorithm 2. In practice, gradient-based optimisers such as L-BFGS-B (Zhu et al., 1997) may settle at local maxima. Thus, we perform optimisation with multiple restarts, with each restart taking a different initial value. A set of variable yet plausible initial values for the restarts can be generated by rerunning Algorithm 2 multiple times with different sliding window widths.

For the GPs defined by eqs. (14), we used squared exponential kernels with constant mean functions (Heinonen et al., 2016). With this assumption, there are six Gaussian process hyperparameters for the model (for each of $L$ and $\sigma$, a mean $\mu$, noise level $\alpha$, and length scale $\beta$). Prior knowledge or a grid search can be used to set these values, although existing work suggests only $\beta$ has a substantial effect on posteriors in most cases (Heinonen et al., 2016). We set $\mu=0$ and $\alpha=1$ (Heinonen et al., 2016), and propose the following method to set $\beta$, which controls how the Gaussian process can change over time. This behaviour is crucial to the adaptivity of the method. If the scale, $\beta$, is too short, the GP will overfit local fluctuations; too large and it will fail to account for real changes in the process over time.

To set the length scale, we used a heuristic based on the expected rate of change of the noise process. Given evenly spaced time points with spacing $\Delta t$ and a user-specified number of time points $N_{c}$, we set $\beta$ as the solution of

for some small value $\zeta=0.01$. This equation imposes that the prior covariance between two values of the Gaussian process $N_{c}$ time points apart is close to zero, thus summarising the prior belief that the noise structure can change over that time scale. Good choices for $N_{c}$ will generally be problem specific. For non-uniform spacing, $N_{c}\Delta t$ could replaced by an appropriate time interval.

In this example, we use the two-parameter logistic growth model:

We demonstrate the results of fitting the non-stationary kernel to synthetic data, generated from a logistic growth model with $r=0.08$, $K=50$, and $f(t=0)=2$ with multiplicative Gaussian noise:

where $y_{i}$ is an observed data point, $f(t;\theta)$ is the ODE model solution, and ${v_{i}\overset{\text{IID}}{\sim}\mathcal{N}(0,\sigma^{2})}$. We set $\eta=2$ and $\sigma=0.0075$. Eqs. (16)&(17) were used to generate eight replicate time series, each with 250 time points. We considered parameter inference for each set of series under three different noise processes: multiplicative (i.e. the true noise process), the non-stationary Laplacian kernel, and IID Gaussian. In each case, Algorithm 1 was used to generate posterior samples from $r$ and $K$. MCMC sampling for model parameters was performed using Pints inference software (Clerx et al., 2019) with three Markov chains and a total of 20,000 iterations on each. The posterior was sufficiently simple here that a non-gradient-based sampler – a form of adaptive covariance algorithm, called the Haario Bardenet method (Haario et al., 2001; Johnstone et al., 2016) – was used opposed to Hamiltonian Monte Carlo. On a desktop processor, each chain took approximately 20 minutes to run. The first half of each chain was discarded as warm-up, and convergence was assessed using the Gelman $\hat{R}$ statistic, requiring $\hat{R}<1.05$ for all parameters (Gelman et al., 2013). To set the GP hyperparameter $\beta$, we used eq. (15) with $N_{c}=200$. The results are shown in Figure 2. In panel (a), the data (from the first replicate) is shown in the top panel, along with the fitted model trajectory. Below, the standard deviation and lag 1 autocorrelation are shown based on the MAP estimates for each replicate and indicate good correspondence with the ground truth. In panel (b), the posterior distributions for the model parameters are shown. The growth parameter, $r$, was most affected by incorrectly assuming IID Gaussian noise, where the IID noise model resulted in estimates with overly inflated uncertainty. This is because model output is most sensitive to $r$ in the first half of the series, where the IID noise model overestimates the noise level. In all cases, the GP method provided a higher fidelity estimate of uncertainty than IID noise; in most cases the location of the posterior is also improved. Another example of the GPs fitted to synthetic data is given in Figure S2. In this example, the true data generating process consists of discrete blocks of different noise models, and the results show the ability of the non-stationary kernel method to find an appropriate smooth approximation.

Instead of assuming that kernel parameters vary smoothly and continuously over time, we now introduce a model that divides the time domain into discrete sections; each with a distinct noise model. A model of this type results in a block structure for the covariance matrix. Assuming the time series is divided into $M$ regimes indexed by $m=1,\dots,M$, and each regime $m$ has a covariance kernel $k_{m}$, the covariance matrix takes the form,

where each block $\Sigma^{(m)}$ is a positive definite matrix, with elements given by ${\Sigma^{(m)}_{ij}=k_{m}(t_{i},t_{j})}$. The covariance matrix from eq. (18) is guaranteed to be positive definite because each block is positive definite. Eq. (18) requires that the autocorrelation between consecutive points at the boundary between two blocks is zero. When autocorrelation varies more gradually, this may be a disadvantage, and could likely be better handled by the non-stationary kernel method in §3. But, for more rapid changes in covariance, the block method should outperform.

To infer the locations of block boundaries is a type of “change point” problem (see, for example, Aminikhanghahi and Cook (2017)), which is a wide field with many methods and fitting processes. We choose a nonparametric approach to change point detection, which has the benefit that the number of boundaries need not be fixed beforehand. In particular, we specify the “restricted product partition model” (PPM) induced by the Pitman-Yor process as a prior on the partitions (Martínez and Mena, 2014). The PPM approach is convenient for learning model parameters, as the posterior over partitions can be inferred jointly with the posterior over ODE model parameters.

To specify the change point model, we place a prior over the partitions which lie in the class of set compositions of the index set $[N]=\{1,2,\dots,N\}$, i.e., $\mathcal{C}_{[N]}$. This class contains all partitions $\{S_{1},S_{2},\dots,S_{M}\}$ such that $S_{m}=\{s_{j}+i:i=1,2,\dots n_{m}\}$, where $s_{0}=0$, $s_{j}=n_{1}+n_{2}+\dots+n_{j}$, and $n_{m}$ is the number of indices in $S_{m}$  (Martínez and Mena, 2014). As an illustrative example, consider the case where $N=3$, in which $\mathcal{C}_{[3]}$ contains the following elements:

$\mathcal{C}_{[N]}$ contains all admissible assignments of the time points to consecutive blocks.

Let $\rho_{N}$ indicate a random variable taking values in $\mathcal{C}_{[N]}$. For a problem with $N$ time points and $k$ regimes, the PPM prior on $\rho_{N}$ with discount parameter $s$ and strength parameter $\phi$ is given by:

In this equation, $x_{n\uparrow}$ denotes the rising factorial or Pochhammer function; i.e., $x_{n\uparrow}=(x)(x+1)\dots(x+n-1)$. Some visual analysis of the prior $p(\rho_{n})$ is presented in Figure S3; $\phi$ serves as a location parameter controlling the number of blocks, while $s$ is a concentration parameter controlling the variance in this number.

We assume a parametric form for all block kernels $k_{m}$, such as the Laplacian kernel given in eq. (8) (which takes two parameters). These parameters are constant within a block, but vary from block to block. The appropriate priors for the kernel parameters will depend on the form of the kernel. For the Laplacian kernel, we place diffuse normal priors on the logarithms of the parameters:

MCMC inference for the posterior is possible using Metropolis-Hastings steps in a split-merge-shuffle sampler, as used in Algorithm 2 of Martínez and Mena (2014). Each iteration of this approach consists of a series of proposals that may be accepted or rejected: first, there is a proposal that affects the number of blocks; this can be either a merge of two consecutive blocks or a split at some random point within an existing block; second, a proposal that keeps the number of blocks the same is used – in this “shuffle” proposal, the boundary between two blocks is randomly moved somewhere within their limits.

In prior work on inference for PPM models (Martínez and Mena, 2014), the MCMC algorithm is restricted to the case where all within-block parameters can be integrated out, allowing calculation of the marginal posterior probability $p(\rho_{N}|y)\propto p(y|\rho_{N})p(\rho_{N})$ for any valid proposed partition $\rho_{N}$. In this work, we relax this assumption, so that we can only calculate the likelihood given explicit values for both the partition and the block parameters. This presents a challenge to inference, as evaluations of the marginal likelihood are essential to calculate the acceptance ratio of the merge and split proposals. In this work, we modify the split-merge-shuffle sampler to handle the non-conjugate case using the ideas of the Sequentially-Allocated Merge-Split (SAMS) Sampler for Dirichlet processes (Dahl, 2005). Briefly, when proposing a split, we simultaneously propose new values for the block parameters in the newly split block using a random walk. When proposing a merge, the current parameter values from the first of the two merged blocks are selected for use in the proposal. At each MCMC iteration, an additional Metropolis-Hastings step is performed to update the explicit block parameter values. The details of our MCMC algorithm are given in §S4.

A key benefit of a nonparametric approach for covariance matrix estimation is that it can learn an appropriate number of blocks from data. This flexibility has a cost, however – in theory, each individual time point could be assigned to its own block. To avoid this outcome, informative priors on $s$ and $\phi$ can be specified. For $s$, we specify the prior: $s\sim\text{beta}(1,1)$. For $\phi$, we place a shifted gamma prior²²2If $X-c\sim\text{gamma}(a,b)$, then $X\sim\text{shifted-gamma}(a,b,c)$. with parameters $(a,b,s)$. We choose the hyperparameters $a,b$ such that, at the prior mean of $s$, prior mass in the number of blocks is heavily concentrated at one block. Appropriate values for $a$ and $b$ may thus depend on the length of the time series and the desired strength of the prior preference for one block; we achieve a moderately informative prior using $a=0.01$ and $b=100$. This specification represents a weak null belief that the noise process is stationary, but allows the noise process flexibility to change if necessary.

We tested the block noise process on synthetic data generated according to the logistic model, eq. (16), with $r=0.08$, $K=50$, and $f(t=0)=2$. In addition, we contrived a noise process with 5 regimes, each of 100 consecutive time points. The first, third and fifth regimes had IID Gaussian noise with $\sigma=3$, the second regime had AR(1) noise with $\rho=0.85$ and $\sigma=3$, and the fourth regime had IID Gaussian noise with $\sigma=30$. A total of 20,000 MCMC iterations were performed using our version of the split-merge-shuffle sampler (see §S4), with the first 10,000 discarded as warm up, and convergence was assessed using the Gelman $\hat{R}$ statistic, requiring $\hat{R}<1.05$ for all model parameters. Each MCMC iteration consisted of one split-merge-shuffle step for the blocks as well as one adaptive covariance step for the model parameters. On this dataset, all iterations took approximately 100 minutes to run on a desktop processor. In the first panel of Figure 3 (a), we show the data and the posterior median of the learned model trajectory, as well as the posterior median of two standard deviations of the inferred noise process about the model trajectory. The next two panels compare the estimated error standard deviation and lag 1 autocorrelation with their true values. This shows that the change point flexible noise process readily captures the five different regimes and learns their boundaries. In Figure 3 (b), we plot the posterior distributions for the logistic growth parameters with three different assumptions for the noise process. The noise process labelled “Correct” indicates an assumption of a multivariate normal likelihood with a block covariance matrix as given in eq. (18), with the block locations and sizes fixed to their correct values. However, in the “Correct” comparator method, the kernel parameters within each fixed block are not known, and are inferred jointly with the ODE model parameters. Meanwhile, the results from the nonparametric block covariance method are labelled with “Block”. For both $r$ and $K$, the IID noise model results in inflated estimates of posterior variance, while both the block covariance method and the true noise model recover precise posteriors.

The human Ether-à-go-go-Related Gene (hERG) encodes the alpha subunit of the potassium channel Kv11.1 that conducts the rapid delayed rectifier potassium current $I_{\text{Kr}}$. This current is of great importance in cardiac electrophysiology and safety pharmacology; reduction of $I_{\text{Kr}}$ by pharmaceutical compounds or mutations can induce fatal disturbances in cardiac rhythm. Interest in this model generally centres on understanding the current response of the hERG channel when a voltage stimuli $V$ is applied. The current can be described with a Hodgkin & Huxley-style structure model (Hodgkin and Huxley, 1952) given by:

where $g_{\text{Kr}}$ is the maximal conductance, and $E_{\text{K}}$ is the reversal potential (Nernst potential) for potassium ions which can be calculated directly from potassium concentrations using the Nernst equation.

The kinetic terms of the model, $a$ and $r$, are governed by:

The model has 9 parameters $\theta=(g_{\mathrm{Kr}},p_{1},p_{2},\dots,p_{8})$ to be inferred, all of which are positive. These parameters are the maximal conductance $g_{\mathrm{Kr}}$ [pS] and kinetic parameters $p_{1},p_{2},p_{3},\cdots,p_{8}$ [s^-1, V^-1, s^-1, $\cdots$, V^-1].

Experimental data of the current are taken from a freely available dataset (Lei et al., 2019a; 2019b), where the voltage stimuli $V$ were designed for parametrising the model.

The logarithm-transformation was applied to all model parameters $\theta$, such that the transformed parameters $\phi=\log(\theta)$ are unconstrained. To account for the impact of this non-linear transformation on the posterior, a Jacobian transformation was applied. Priors for $\phi$ were selected using existing literature results (Lei et al., 2019a; 2019b), and, for each element of $\phi$, a weakly informative prior Gaussian distribution was used (see Table S1 for the prior hyperparameters).

For six different cells, the model parameter posteriors were obtained via MCMC using the IID noise model and the non-stationary Laplacian kernel flexible noise model. To obtain posterior samples, the simulated tempering population MCMC algorithm was used (Jasra et al., 2007), with convergence assessed using the Gelman $\hat{R}$ statistic (Gelman et al., 2013), and the first half of each chain discarded as warm up. For the non-stationary Laplacian kernel, we used Algorithm 1 for inference and Algorithm 2 for initialisation. The fits for each of six cells are shown with the time series data in Figure S4.

Figure 4 shows the central 95% posterior distribution ranges for all nine model parameters, assuming either IID Gaussian noise (horizontal axis) or the non-stationary noise process (vertical axis). There were significant differences in the parameter estimates for almost all parameters, with much of probability mass not overlapping the IID=Laplacian line. Additionally, the more sophisticated noise model resulted in substantially higher posterior variance for several model parameters, notably including $g_{\text{Kr}}$, $p_{2}$ and $p_{4}$. Cell A04 is an outlier: this is likely because this cell has a region of drastic misspecification in much of the time series, from $t=6$ to $t=10$. While the model fits for all six cells indicate short regions of misspecification, which is particularly apparent after the drops in current around $t=2$ and $t=14$, cell A04 (and to a lesser extent, A07) suffer from more extensive misspecification. The data and inferred fits for cell A04 are shown in Figure 5. The non-stationary noise model detects the central misspecified region by assigning high variance and autocorrelation in the middle of the time series. In the time series for cells A04 and A07, the poor fit between model and data may be largely explained by the fact that our model in this study (§6.1) fails to account for experimental artefacts in the voltage clamp experiment, such as leakage current—these artefacts may explain much of the cell-to-cell variability observed in these experiments (Lei et al., 2020). Thus, the high levels of standard deviation and autocorrelation detected in these time series suggest that a more detailed model of the experiment is necessary in order to understand these cells and correct the regions of obvious poor fit.