Bayesian Inference for Multiple Datasets

1. Introduction

A fundamental problem in modelling is the estimation of unknown parameters from empirical data . Frequently however, practical situations require estimation of individual parameters for a large number of datasets. Fitting the same model to multiple datasets and extracting parameters specific for each dataset is becoming common when dealing with biological and epidemiological data. Examples range from estimation of specific parameters across different genotypes [1], cell lines [2]; patients [3,4], health districts or country specific parameters [5,6]. The number of individual datasets can increase up to hundreds of thousands when linking geostatistical maps and disease transmission models [7,8,9]; or combining remote sensing data with crop growth models [10,11]. Inference would be prohibitively slow when using methods such as Markov chain Monte Carlo sampling [12] or importance sampling [13]. Recently, a new approach, called adaptive multiple importance sampling (AMIS), has shown promising results when dealing with sampling from multiple targets simultaneously [8,14,15].

In this work we use the adaptive multiple importance sampling approach to fit von Misses distribution to wind trajectories. The von Mises distribution is commonly used for modelling circular data [16], as well as for statistical modelling of wind direction [17]. The von Mises distribution has two parameters: the directional mean and concentration. Temporal changes in the directional mean can show seasonal variation in global wind patterns, while the concentration can be used as a measure of stochasticity associated with natural variation in turbulent wind flow. Values of concentration close to zero indicate a uniform angular distribution, and high values imply low variance that is consistent with stronger directional flow. Knowing wind direction and the uncertainties related to it are important when evaluating risk or planning mitigation actions associated with movements of pest and disease. For example, visualization of long-range atmospheric transport have been utilized to improve risk estimates and early warning in operational crop disease and insect pest management systems [18,19].

The Lagrangian atmospheric dispersion model NAME can forecast or reconstruct wind trajectories at a high spatial and temporal resolution[20]. Originally developed by the UK Met Office to model the release, transport, dispersion and removal of hazardous material in the atmosphere [20], NAME has subsequently been adapted to model long-distance wind-borne transport of insect vectors of viral disease [21] and fungal spores [18]. Furthermore, the Lagrangian atmospheric dispersion model NAME can simulate wind trajectories that are stochastic so that an ensemble of trajectories provides a probabilistic distribution of the dispersal of particulates such as spores or insects [20]. For our analysis we considered a domain extending across five sub-Saharan countries (Kenya, Ethiopia, Somalia, Eritrea and Djibouti). This area was selected because it was affected by the 2019-2021 desert locust upsurge [22]. First, to gauge the computational effort required, we use importance sampling to recover parameter values from synthetic data for a single dataset. Then using synthetic data and three datasets we show that AMIS can significantly reduce computational effort while maintaining the quality of inference. Finally, we test the performance of the method on a very large number of datasets (>11,000) using wind trajectories derived from a Lagrangian Particle Dispersion Model driven from a 3D meteorological data.

2. Materials and Methods

2.1. The von Mises Distribution

A random variable $\theta \in [0,2\pi ]$ has a von Mises distribution, if its probability density function is defined as [23]:

$$\begin{array}{c}\hfill f_{VM}(\theta ,\mu .\kappa )=\frac{1}{2\pi I_0\left(\kappa \right)}exp\left[\kappa cos\left(\theta -\mu \right)\right].\end{array}$$

(1)

Here $I_0\left(\kappa \right)$ is the modified Bessel function of the first kind and order zero, $\mu \in [0,2\pi ]$ is the directional mean, and $\kappa \ge 0$ is the directional concentration.

2.2. Synthetic Wind Data Used for Initial Test

We used function rvonmises from the R package Rfast [24] to sample random variables from the parameterised von Misses distribution.

2.3. NAME Wind Trajectories Based on Historic Meteorology Data

The Lagrangian atmospheric dispersion model NAME was used to reconstruct wind trajectories at a high spatial ( 10 km) and temporal (3h) resolution [20], relative to the domain of interest extending across multiple countries. These NAME simulations are driven using historic analysis meteorology data from the global configuration of the Unified Model (UM), the Met Office’s operational Numerical Weather Prediction model [25]. NAME wind trajectories were calculated starting from 3860 source locations spaced on a regular 20km × 20km grid across Ethiopia, Somalia, Eritrea and Kenya (Figure 1(A)). For a given date and location source, 1000 individual trajectories were calculated, each commencing 2h after local sunrise and terminating 1h before local sunset (Figure 1(B)). There is stochastic variation in the direction of individual realisations starting from the same location, due to the turbulent nature of the atmosphere as represented in NAME. We calculated the angle for each trajectory between a vector corresponding to the East direction and a vector connecting the start and end location of a trajectory for a specified time.

2.4. Bayesian Inference

2.4.1. Importance Sampling (IS)

Importance sampling is based on using weighted samples from a proposal distribution [26]. The corresponding importance weights for a sampled set of parameters $({\mu }_i,{\kappa }_i)$ can be calculated as:

$$w_j^i=\frac{f_{VM}\left({\theta }_j,{\mu }_i,{\kappa }_i\right)}{\varphi ({\theta }_j,{\mu }_i,{\kappa }_i)}$$

(2)

where $f_{VM}\left({\theta }_j,{\mu }_i,{\kappa }_i\right)$ is the target density function (i.e. the von Mises distribution function) and $\varphi (\theta ,\mu ,\kappa )$ is the proposal density function.

2.4.2. Effective Sample Size (ESS)

We use Kish’s effective sample size [27] as a measure of the quality of the posterior distribution. We denote the ESS for dataset i after iteration l as,

$$\begin{array}{c}\hfill ES{S}_i={\left(\sum _{j=1}^N{\left({\widehat{w}}_i^t\right)}^2\right)}^{-1}.\end{array}$$

(3)

Here ${\widehat{w}}_i$ is a normalised weight, so that ${\sum }_i^n{\widehat{w}}_i=1$.

2.4.3. Adaptive Multiple Importance Sampling (AMIS)

The Adaptive Multiple Importance Sampling algorithm (AMIS) is a technique that recycles samples from all previous iterations when constructing a proposal distribution or calculating weights [28]. For each iteration $l=1,\dots ,L$, we sample $N_l$ parameter sets. The corresponding importance weights for a sampled set of parameters $({\mu }_i,{\kappa }_i)$ can be calculated as:

$$\begin{array}{c}\hfill w_j^i=\frac{f_{VM}\left({\theta }_j,{\mu }_i,{\kappa }_i\right)}{\frac{1}{N_1+\dots N_L}{\sum }_{l=1}^L{N}_l{\varphi }_l({\theta }_j,{\mu }_i,{\kappa }_i)}\end{array}$$

(4)

where ${\varphi }_l(\mu ,\kappa )$ is a proposal distribution at iteration l.

For the AMIS algorithm, we set a required ESS for all of the pixels, denoted by $ES{S}^R$. We call datasets ’active’ if they have an ESS below the target, and use the weights for the active datasets to design the proposal for the next iteration of the AMIS algorithm. This targets sampling towards areas of the parameter space that benefit datasets that have not yet reached the required ESS. At iterations $l=2,...$, we use the mean weight over the active dataset to determine the next proposal distribution. A suitable proposal can be found by fitting a mixture of Student’s t-distributions to the weighted samples $({\theta }_j,{\overline{w}}_j^i)$, for $j\in \{1,\cdots N_1+\cdots +N_i\}$. The algorithm continues until all datasets meet the ESS requirement or the maximum number of iterations, I, has been completed.

3. Results

3.1. Illustrative Simulations

3.1.1. Parameter Estimation Using Importance Sampling for a Single Dataset

First, we sampled 1000 angles using parameter values $\mu =\pi$ and $\kappa =5$ (Figure 2(A)). For the IS, we chose a uniform distribution function as our proposal: $\varphi (\theta ,\mu ,\kappa )\sim U[0,360]\times U[0,700]$. We varied the number of parameters sets sampled from the proposal from ${10}^3$ to $2\times {10}^6$ and calculated the effective sample size (ESS) as a function of the number of parameter sets sampled Importance sampling produced a reasonable posterior distribution that recovered the original parameter values (Figure 2(B)-(C)) and provided a corresponding fit of the von Mises distribution function to the original (Figure 2(D)). Although the ESS increased with sample size, the ESS was >25 for $2\times {10}^6$ samples for this single dataset (Figure 2(E)).

3.1.2. Parameter Estimation Using AMIS for Multiple Datasets

Next, we considered three datasets with differing means and variances. We sampled 1000 angles with the following values of parameters: set 1 $({\mu }_1,{\kappa }_1)=(0.52,100)$; set 2 $({\mu }_2,{\kappa }_2)=(1.9,20)$; and set 3 $({\mu }_3,{\kappa }_3)=(4.2,5)$ (Figure 3(A)). We chose these parameters to produce angle distributions which do not overlap, as well as having different dispersion around their mean. For this case study, we have fitted parameters for the three sets simultaneously using the AMIS algorithm. As in the previous example, we chose a uniform distribution function as our proposal: $\varphi (\theta ,\mu ,\kappa )\sim U[0,360]\times U[0,700]$. We set the required effective sample size $ES{S}_R=100$ and we sampled 1000 parameter sets each iteration. The algorithm stopped when all datasets had $ES{S}_i\ge ES{S}_R$. The first dataset reached the required ESS after 6000 sampled parameter sets (Figure 3(B)). The other two datasets reached the required ESS after 8000 parameter sets to achieve the required $ES{S}_R$. The posterior distributions recovered the original parameters (Figure 3(C)).

3.2. Application to NAME Wind Trajectories

Finally, we simultaneously fitted the von Mises distributions to the NAME wind trajectories. We chose three dates: 15 January 2020, 15 August 2020 and 15 January 2021. This resulted in the estimation of parameters for 11,580 datasets, corresponding to 3 x 3860 datasets on the 20km × 20km grid. We chose a uniform distribution function as our proposal: $\varphi (\theta ,\mu ,\kappa )\sim U[0,360]\times U[0,700]$. We set the required effective sample size $ES{S}_R=100$ and we sampled 1000 parameter sets each iteration. It took about $2*{10}^6$ sampled parameters for all datasets to reach the required ESS (Figure 4(A)). For each grid cell and the three days we provide 100 parameter sets sampled from the posterior distributions in the supplementary materials.

To illustrate global wind trends, we have drawn a single parameter sample from the posterior distribution for each of the 3860 sites and plotted the spatial distribution of these parameters (Figure 4(B)) for each of the three selected dates. The directional mean for January had similar trends for both 2020 and 2021, i.e. most of wind trajectories had southwards orientations, but some differences in direction could be seen for small areas. The concentration values showed more variation, with northerly wind concentrated less in 2021 than in 2020. Meanwhile, wind trajectories were directed towards the north on 15 August 2020, with a larger range of concentrations in comparison with wind trajectories in January for both 2020 and 2021.

4. Discussion

In situations when the number of datasets is large, performing Bayesian inference can be time consuming. Here we have investigated how adaptive multiple importance sampling can be used to effectively estimate parameters for multiple datasets. As an illustration of the methodology, we fitted the von Mises distribution to wind direction data. For a comparison, we used importance sampling to fit parameters to a single wind angle distribution using synthetic data. Inference using importance sampling recovered the parameter values and captured the original distribution for wind trajectories, but it was prohibitively slow and computationally inefficient ($2\times {10}^6$ parameter sets to achieve ESS=25). We showed that AMIS significantly reduced computational effort while maintaining the quality of inference (8000 parameter sets to achieve ESS=100) for three independent datasets. Finally, our results showed that AMIS was able to achieve ESS=100 for >11,000 location and date specific parameters for the NAME wind trajectories with similar computational requirements as a single dataset fitted using importance sampling ($2\times {10}^6$ parameter sets).

The migration pathways achieved by flying insects are highly dependent on the seasonal variation of wind patterns [29,30,31]. For example, swarms follow seasonal shifts of the Intertropical Convergence Zone (ITCZ), which can take them to areas of rainfall [32,33,34,35,36]. We fitted the von Mises distribution to NAME wind trajectories at two different seasons (January and August). Our analysis confirmed the influence of the ITCZ on local wind trajectories, i.e. southwards orientation in January and northwards orientation in August. Certain areas tended to have a smaller concentration in the direction of wind trajectories, and hence greater variability which is important to account for when predicting long-distance windborne dispersal. Another application of fitting the von Mises distribution to NAME wind data is to compare annual differences in wind behaviour at high resolution, as illustrated in Figure 4(B).

As the amount of data continues to increase at an exponential rate [37] and more multi-set data become available, it will be crucial to have computationally efficient methods for parameter estimation. Our proposed methodology can be easily applied to any type of a problem where the target distribution can be formulated due to the universal way the proposal density function is constructed at each iteration. Further research is required to guide the optimal choice of ESS and number of sampled parameters at each iteration for a particular type and size of datasets.