# Inference and model determination for Temperature-Driven non-linear
Ecological Models

Marios Kondakis  Nikolaos Demiris and Ioannis Ntzoufras  Athens University of
Economics and Business, Athens, Greece.  and Nikos E. Papanikolaou
Agricultural University of Athens, Athens, Greece.  Greek Ministry of Rural
Development and Food, Athens, Greece.

###### Abstract

This paper is concerned with a contemporary Bayesian approach to the effect of
temperature on developmental rates. We develop statistical methods using
recent computational tools to model four commonly used ecological non-linear
mathematical curves that describe arthropods’ developmental rates. Such models
address the effect of temperature fluctuations on the developmental rate of
arthropods. In addition to the widely used Gaussian distributional assumption,
we also explore Inverse Gamma–based alternatives, which naturally accommodate
adaptive variance fluctuation with temperature. Moreover, to overcome the
associated parameter indeterminacy in the case of no development, we suggest
the Zero Inflated Inverse Gamma model. The ecological models are compared
graphically via posterior predictive plots and quantitatively via Marginal
likelihood estimates and Information criteria values. Inference is performed
using the Stan software and we investigate the statistical and computational
efficiency of its Hamiltonian Monte Carlo and Variational Inference methods.
We explore model uncertainty and use Bayesian Model Averaging framework for
robust estimation of the key ecological parameters.
_keywords:_ Arthropods’ Developmental rates, Bayesian Model Averaging,
Ecological models, Marginal likelihood estimation, Model comparison, Bayesian
computation

00footnotetext: Address of correspondence: Marios Kondakis, Department of
Statistics, Athens University of Economics and Business, 76 Patission str., PO
10434, Athens, Greece.
<EMAIL_ADDRESS>
##  1 Introduction

Studying the population evolution of arthropod pests, as well as of biological
control agents, is of great importance for the crop primary production,
agricultural infrastructures, spreading diseases and consequently the economy
(Bradshaw et al., 2016). Temperature and body size are two major determinants
that influence the metabolic, survival, growth and reproduction rates which
control the ecological processes at all levels of arthropods’ life (Brown et
al., 2004). Biological control is facilitated when the climate responses of
biocontrol agents are understood, especially to temperature. The thermal
thresholds for insect development can be estimated using several functional
forms (Kontodimas et al., 2004).

In all but the simplest cases, mathematical modelling is an indispensable tool
for understanding the resulting developmental scheme (Kontodimas et al.,
2004). However, fitting ecological models for developmental rates is not
straightforward, typically because the mathematical forms are not linear
(Papanikolaou et al., 2019) and the actual biochemical reactions of insects or
environmental factors responsible for their growth may remain unobservable.
Therefore, we adopt the Bayesian paradigm to population dynamics’ modeling and
inference since it naturally accounts for latent parameters and their
uncertainties. Nonetheless, there are significant challenges in designing
statistical methods that work efficiently in a wide range of ecological
applications. Stan (Carpenter et al., 2017) provides a BUGS-like interface to
model building and the ability to run it via different languages and operating
systems.

###  1.1 Motivation for modeling developmental rate

Insects and mites, as ectotherms, regulate their body temperature according
the environment they live in (Norris and Kunz, 2012). This affects the rate of
metabolism, i.e., the biochemical reactions that allow the processes of
production and release of energy, as well as the synthesis of necessary
molecules that serve as structural or functional components (Neven, 2000). In
fact, temperature affects the functionality of enzymes, which act as catalysts
for these biochemical reactions (van der Have, 2008). Consequently, within a
range of temperatures in which insects and mites develop and reproduce,
various biological features are affected (Broufas and Koveos, 2001; Huey and
Berrigan, 2001; Broufas et al., 2007; Nedvěd, 2009; Papanikolaou et al., 2013,
2014). Thus, their performance is indebted to several temporal fluctuations in
terms of population size through time. The empirical finding of the initial
increase in the growth rate of insects and mites in relation to temperature,
followed by its sharp decline, formed the basis for the development of various
mathematical models of its description (Kontodimas et al., 2004). These models
allow the estimation of the lower and upper thermal limits, i.e. the lowest
and highest temperature, respectively, at which the growth rate is zero, as
well as the temperature at which it receives its maximum value. Understanding
populations’ growth rate is of importance, as their assessment can lead to
decisions on their management (Hare et al., 2011), particularly under the
pressure of climatic change (Bradshaw et al., 2016).

###  1.2 Historical overview with models and techniques used in literature

This paper investigates some popular non-linear ecological models that
describe the rate of insects’ and mites’ development within a Bayesian
context. We explore the computational and statistical efficiency of
Hamiltonian Monte Carlo (HMC) (Betancourt, 2017; Neal, 2010) and Variational
Bayes (VB) (Blei et al., 2017), a challenging task in the present setting due
to distinct features of the entertained models, including truncation. Both the
widely used Gaussian distributional assumption and a newly-developed Inverse
Gamma-based version are explored in order to model the developmental rate
distribution of insects and mites. We compare the models using information
criteria, models marginal likelihood estimates and graphical tools. In
addition, model averaging techniques are used to provide robust estimates of
the parameters of interest. A distinct feature of these models is that zero
count data make one parameter of interest indeterminable and model fit
potentially misleading. We propose ways to overcome this indeterminacy by
applying the Zero Inflated Inverse Gamma distribution while carefully
connecting the probability of non-zero development to the predictor.

The remainder of this paper is structured as follows. The next section
contains the ecological models we develop. Section 3 investigates two separate
real-life examples without and with zero-rates, and Section 4 reports our
findings, while the paper concludes with discussion.

##  2 Ecological models

A variety of different linear and non-linear equations have been used to
describe the rate of development of insects and mites and to estimate their
thermal limits. Such linear approximations enable the calculation of lower
developmental threshold and thermal constant within a small temperature range,
usually 15-30oC (Campbell et al., 1974; Wagner et al., 1991; Jarošík et al.,
2002; Kontodimas et al., 2004). However, the relationship between development
and temperature becomes non-linear outside that range. Thus, in order to
accurately predict developmental rates across the spectrum, the use of non-
linear ecological models is required (Wagner et al., 1991; Kontodimas et al.,
2004; Damos and Savopoulou-Soultani, 2012).

###  2.1 Non-linear ecological models

In the present work, four commonly used non-linear ecological models are
considered. Specifically, the Bieri (Bieri et al., 1983), the Briere (Briere
and Pracros, 1998; Briere et al., 1998, 1999), the Analytis (Analytis, 1981)
and the Lactin (Lactin et al., 1995) models are implemented. Developmental
time is the duration between life stages of the insects or mites (Wagner et
al., 1991). The response variable $y(T;{\theta})$ describes the developmental
rate and it is defined as the reciprocal of the days until the completeness of
a particular developmental event. Herein, $T$ denotes the predictor variable,
the absolute temperature measured in Celsius degrees, while ${\theta}$ denotes
the parameter vector of the model. Typically, the expected developmental rate,
$r(T;{\theta})$ is modelled and the four aforementioned models are presented
below.

####  2.1.1 Bieri Model

In the Bieri model, the developmental rate is defined as

$r(T;{\theta})=\alpha\cdot(T-T_{m_{1}})-{\beta}^{(T-T_{m_{2}})}$ (1)

where $\alpha$, $\beta$, $T_{m_{1}}$ and $T_{m_{2}}$ are the model parameters.
In particular, the values of $T_{m_{1}}$ and $T_{m_{2}}$ lie close to the real
lower and upper thermal thresholds $T_{min}$ and $T_{max}$, at which the
development starts or ceases respectively. The exact values of the thermal
thresholds are derived implicitly as the lower and upper roots respectively of
the response function. The parameter $\alpha$ corresponds approximately to the
rate of increase in the linear model at vital temperatures (Bieri et al.,
1983). The response variable in (1) is a concave function of the temperature.
Parameter $\alpha$ is defined in the interval $\left(0,1\right)$ while $\beta$
determines the decrease of the developmental rate at higher temperatures
(Bieri et al., 1983) when $\beta$ exceeds unity. According to (1), we observe
that:
$r(T_{m_{1}})=-{\beta}^{(T_{m_{1}}-T_{m_{2}})}$ and
$r(T_{m_{2}})={\alpha}(T_{m_{2}}-T_{m_{1}})-1$,
which suggests that $r(T_{m_{1}}),r(T_{m_{2}})\in(-1,0)$ and leads to the
following inequality

$r(T_{m_{1}})<r(T_{min})<r(T_{max})<r(T_{m_{2}}).$

The temperature at which the maximum developmental rate occurs is called
optimum and it is denoted by $T_{opt}$. In the Bieri model, it is given by

$T_{opt}=T_{max}+\frac{\log{\alpha}-\log\left({\log{\beta}}\right)}{\log{\beta}}.$

In Fig. 1 some curves are generated by the Bieri model in (1) when
$T_{m_{1}}=5^{\circ}$C and $T_{m_{2}}=35^{\circ}$C while the other parameters
vary.

Figure 1: Bieri developmental rates

####  2.1.2 Briere Model

The Briere model is the most popular and parsimonious model. The developmental
rate is defined as

$r(T;{\theta})=\begin{cases}\alpha\cdot
T\cdot(T-T_{min})\cdot\sqrt{(T_{max}-T)}&\mbox{ for }\,T_{min}<T<T_{max}\\\
0&\mbox{ otherwise }\\\ \end{cases}$ (2)

where $\alpha$, $T_{min}$ and $T_{max}$ are model parameters. Particularly,
$T_{min}$ and $T_{max}$ are exactly the lower and upper thermal thresholds at
which the development starts or ceases respectively while parameter $\alpha$
is an empirical constant (Briere et al., 1998). The response variable in (2)
is again a concave function of the temperature. Parameter $\alpha$ is defined
in $\left(0,1\right)$ whereas the existence of the square root in (2) ensures
that the developmental rate declines sharply at higher temperatures. The
optimum temperature in the Briere model is given by

$T_{opt}=\frac{1}{10}\left\\{4\cdot T_{max}+3\cdot T_{min}+\sqrt{{(4\cdot
T_{max}+3\cdot T_{min})}^{2}-40\cdot T_{min}\cdot T_{max}}\right\\}.$

In Fig. 2 some curves are created by the Briere model in (2) when
$T_{min}=5^{\circ}C$ and $T_{max}=35^{\circ}C$ while the parameter $\alpha$
varies.

Figure 2: Briere developmental rates

####  2.1.3 Analytis Model

The developmental rate in Analytis model is defined by

$r(T;{\theta})=\begin{cases}\alpha\cdot{(T-T_{min})}^{n}\cdot{(T_{max}-T)}^{m}&\mbox{
for }\,T_{min}<T<T_{max}\\\ 0&\mbox{ otherwise }\\\ \end{cases}$ (3)

where $\alpha$, $T_{min}$, $T_{max}$, $n$ and $m$ are parameters of this
model. The exponents $n$ and $m$ in (3) are empirical constants that determine
the rate of growth and decrease of the developmental rate respectively
(Analytis, 1981). Both these parameter take values in $(0,+\infty)$ but in
order to reduce the computational burden we may restrict them in a subset of
the form $(0,c)$, for some constant $c>0$. Finally, $\alpha$ takes values in
$\left(0,1\right)$ interval. The optimum temperature in this model is given by

$T_{opt}=\frac{n\cdot T_{max}+m\cdot T_{min}}{n+m}.$

The Analytis model has a multiplicative polynomial structure in which the
exponents change as parameters to be estimated. Such a structure needs some
empirical driven tuning when defining its thermal parameters space.
Especially, in the case of one dataset, we assumed that $T_{min}$ is greater
than ${4^{\circ}C}$. Some curves generated by the Analytis model in (3) are
depicted in Fig. 3 when $T_{min}=5^{\circ}C$ and $T_{max}=35^{\circ}C$, while
parameters $n$ and $m$ vary.

Figure 3: Analytis developmental rates

####  2.1.4 Lactin Model

The Lactin model includes four parameters and the developmental rate is
defined as

$r(T;{\theta})=\lambda+e^{\rho\cdot T}-e^{\rho\cdot
T_{m}-\frac{(T_{m}-T)}{\varDelta}}$ (4)

where $T_{m}$ is associated with the upper thermal threshold $T_{max}$ since
it tends to this value when $\lambda$ tends to be zero. Parameter $\lambda$
represents an asymptotic level of the developmental rate value in (4) that is
approximated when predictor $T$ tends either to $-\infty$ (extremely low
temperatures) or to the threshold parameter $T_{m}$. Thus, in the event that
the $\lambda$ is non-negative, the Lactin ecological function does not have a
lower thermal threshold ($T_{min}$) and at the same time $T_{m}\leq T_{max}$
as the developmental rate value in (4) is limited above zero level, at which
the maximum thermal thresholds undergoes. In the case that $\lambda$ is
negative, $T_{m}>T_{max}$ and $T_{min}$ sample space is in the
$\left(-\infty,T_{max}\right)$ interval. Parameter $\varDelta$ is positive and
it determines the descent steepness of the developmental rate. It expresses
the temperature range between the value at which the response function begins
to descend and the value of the $T_{m}$ parameter. When $\varDelta$ is less
than one, the rate of descent is very high, although in the other case the
rate of descent is lower and similar to the other ecological models. This
feature triggers discontinuity and lack of fit problems. Hence, in this work
we define $\varDelta$ in $(1,+\infty)$ interval in order to avoid such
problems. Parameter $\rho$ describes the acceleration of the function from low
temperatures to the optimal temperature (Lactin et al., 1995). The response
function in (4) has one inflection point, a maximum point at the optimum
temperature and asymmetry about this point (left skewed). Also the function
has a sharp drop after the optimum temperature, which is achieved by setting
$\rho$ in $(0,{\varDelta}^{-1})$. The actual upper thermal threshold $T_{max}$
is evaluated as the higher root of response function in (4). The optimum
temperature in the Lactin model is given by

$T_{opt}=T_{m}-\frac{\varDelta\cdot\log{(\rho\cdot\varDelta)}}{\rho\cdot\varDelta-1}.$
(5)

In addition, the temperature at the inflection point of the Lactin curve is
given by

$T_{inf}=T_{opt}-\frac{\varDelta\cdot\log{(\rho\cdot\varDelta)}}{\rho\cdot\varDelta-1}.$
(6)

In Fig. 4 some curves are created by the Lactin model in (4) when
$T_{m}=35^{\circ}C$ and the parameters $\Delta$, $\rho$ and $\lambda$ vary.

Figure 4: Lactin developmental rates

###  2.2 Ecological features of the models

There are some basic common features in all of the above-mentioned ecological
models.

There is no growth below the lower temperature threshold $T_{min}$ or above
the upper temperature threshold $T_{max}$. Specifically, in the case of the
Briere and Analytis models, the developmental rate is positive and is defined
only between the two parameters of the thermal thresholds. Also, the
developmental rate in the Bieri and Lactin models can take negative values
that cannot be interpreted. To accommodate these characteristics, we use
initial values of the parameters that allow the ecological function to receive
positive values.

The developmental rate is an asymmetric curve left skewed of its maximum
point. It increases and reaches a maximum at optimal temperature while it
declines rapidly down to zero at the higher temperature threshold $T_{max}$
that is considered as lethal temperature. It includes an inflection point,
with the exception of the Bieri model which is a concave function of the
temperature. The structure of the Lactin model makes it susceptible to a type
of exponential pattern in data-values (Fig. 4). The Briere model, on the
contrary, has a specific structure, which has the first derivative of class
$\mathcal{O}\left(T^{2.5}\right)$ and can hardly trace exponential changes in
data values (Fig. 2). Also, the model of Analytis has polyonimic structure as
the model of Briere but it does contain more degrees of freedom because the
exponents it includes are unknown parameters (Fig. 3). The Bieri model also
adopts an exponential reduction after the thermal optimum threshold but can
only follow the developmental increase of the dataset linearly (Fig. 1). All
four ecological models have been used in the literature to provide reasonable
estimates of the thermal thresholds of several anthropods’ developmental rates
at various stages (Bieri et al., 1983; Kontodimas et al., 2004; Aghdam et al.,
2011). However, the Bieri model is underutilized in the literature, so we
include it in this study to gain a better understanding of ecological models
that describe temperature-dependent development. Additionally, in (Kontodimas
et al., 2004; Aghdam et al., 2011) the ecological models were compared based
on the accuracy of the real data thermal threshold estimates, the adjusted
coefficients of regression ($R^{2}$), and the residual sum of squares values.
In summary, the Briere and Analytis models, appear to overestimate and
underestimate the upper and lower thermal thresholds, respectively, whereas
the Bieri and Lactin models appear to meet the majority of the criteria used
in the comparisons. Despite these minor differences, all of the above models
appear to provide higher $R^{2}$ values than other models in real-data
applications in the literature and we include all four in the current work.

###  2.3 Measurement Error

Probabilistic random error due to chance and systematic error due to data with
excessive zeros is added in order to include uncertainty in the ecological
models already listed. In section 2.3.1, we present the notation adopted and
the probability schemes implemented in this analysis.

####  2.3.1 The likelihood of the data

Let $y_{i}$ represent the observed developmental rates of the $i^{th}$
individual observed at $T_{i}$ temperature, where $i=1,\dots,N$. We consider
$y_{i}$ to be independent response variables counted as the reciprocal of the
number of days until the development of the $i^{th}$ individual takes place
and the range of its value lies in the $(0,1)$ interval. Furthermore, let $y$,
$T$ and $\theta$ be the vectors of the response, the predictor and the
parameters respectively then the conditional expectation $E(y|T,\theta)$ is
considered to vary according to the ecological function $r(T;\theta)$
presented in (1), (2), (3) and (4) respectively.
The data distribution is denoted by

$p({y}|{T,\theta})=\prod_{i=1}^{N}p(y_{i}|T_{i},\theta).$

In this study we consider the Gaussian and the Inverse-Gamma distributions as
the distribution of the response data $y_{i}$. Specifically, the Gaussian
distribution is used broadly in the literature as a good approximation to most
unimodal distributions with finite variance due to the general form of the
Central Limit Theorem. Thus, it can be used as to approximate a more
complicated model likelihood of the data. The non-linear model of the
independent response rates has the Gaussian distribution given by

$y_{i}\sim N(r(T_{i};{\theta)},\sigma^{2})$ (7)

where the mean of the Gaussian likelihood is driven by the respective
ecological model, whereas its standard deviation $\sigma$ is considered as an
unknown parameter. We consider a weakly informative prior distribution for
$\sigma$ like: the Inverse-Gamma
$\text{Inv}\Gamma\left({10}^{-3},{10}^{-3}\right)$. The non-linear model used
in (7) is a constant variance model among different temperatures. In addition,
it allows for zero observed rate values, which occur when the stage of the
insect does not change in perpetuity. Such modeling, however, has the
disadvantage of lack of interpretability in the case of estimated negative
rate values.

On the other hand, the Inverse-Gamma distribution is a plausible alternative
for modeling positive observed rates, as it handles positive values that
describe ratios such as developmental rates whose inverse are positive counts
(such as days passed until the expected development occurs) that can be
described by the Gamma distribution. The non-linear model of the independent
response rates has the Inverse-Gamma distribution given by

$y_{i}\sim\text{Inv}\Gamma(\zeta,(\zeta-1)\cdot r({T_{i}};{\theta}))$ (8)

where $\zeta$ is the shape parameter of the Inverse-Gamma distribution. The
mean of the distribution equals $r({T_{i}};{\theta})$, whereas the variance
equals $\frac{r({T_{i}};{\theta})^{2}}{(\zeta-2)}$. The Inverse-Gamma
likelihood in (8) is a natural alternative to model observed positive rates.
Herein, its mean is driven by the respective ecological function, whereas its
variance depends both to the shape parameter $\zeta$ and the ecological
function $r({T_{i}};{\theta})$ as well, which allows the variance to be
temperature dependent. For the prior distribution of the shape parameter
$\zeta$ we have chosen a weakly informative
$\Gamma\left({10}^{-1},{10}^{-2}\right)$.

####  2.3.2 Zero-rates case

There are situations in which any insect development does not occur throughout
the cohort study. This is indicated by zero values in the response variable,
which can theoretically be interpreted as the number of days required for the
insect to move to its next stage never ending, implying that the developmental
rate is the reciprocal of infinity. In such cases, the MCMC sampling procedure
can be extended either to include prior information about the case of no
development by adjusting the prior knowledge of the parameters concerned or to
include a zero-inflation scheme. For this work, we suggest the use of a Zero
Inflated Inverse Gamma distribution, which gives zero value with probability
$p_{i}$ for the $i^{th}$ insect observed at $T_{i}$ temperature. Particularly,
the probability density function of the observation $y_{i}$ is:

$\begin{split}P(y_{i}|T_{i},\theta)=\left\\{\begin{matrix}p_{i}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\textrm{
if }y_{i}=0\\\ P_{\text{Inv}\Gamma}\left\\{\zeta,(\zeta-1)\cdot
r(T_{i};{\theta})\right\\}\cdot\left(1-p_{i}\right)\textrm{ if }y_{i}\neq
0\end{matrix}\right.\\\
\textrm{logit}\left(1-p_{i}\right)=\left[c\cdot\left\\{r(T_{i};{\theta})-k\right\\}\right]\Rightarrow
p_{i}=\frac{1}{e^{\left\\{c\cdot\left(r(T_{i};{\theta})-k\right)\right\\}}+1}\end{split}$
(9)

where $c$ is a constant positive parameter that has the opposite order of
magnitude of the sample mean $\overline{y}$, while $k$ is the inflation point
of logit link function where the probability of zero $p_{i}$ is equal to
$\frac{1}{2}$. The constant $k$ is associated with the constant $c$ and can be
chosen so that the zero rate of development matches the probability $p_{i}$ at
a predetermined level like 0.9. According to the real data example presented
below in the results section, the proposed values that satisfy the above
criteria for constants $c$ and $k$ are $10^{2}$ and $5\cdot 10^{-3}$,
respectively.

As the developmental rate of $r(T_{i})$ increases to its maximum, the
probability $p_{i}$ in (9) decreases towards zero. On the other hand, when
$r(T_{i})$ tends to be a very small number, the probability $p_{i}$ in (9)
tends to be one.

####  2.3.3 Priors

In the Bieri model (1) the prior distribution of the parameters $\alpha$,
$\beta$, $T_{m_{1}}$ and $T_{m_{2}}$ are shown in Table 1 respectively.

In the Briere model (2), the transformation $\widetilde{a}=-\log(\alpha)$ is
considered instead of the original parameter $\alpha$. Weakly informative
priors are considered for $\widetilde{a}$, $T_{min}$ and $T_{max}$ as shown in
Table 1 respectively. In the Analytis model (3) the transformation
$\widetilde{a}=-\log(\alpha)$ is considered instead of $\alpha$ with the
weakly prior distribution $\Gamma\left({10}^{-1},{10}^{-2}\right)$.
Additionally, the prior distributions utilized for parameters $m$, $n$,
$T_{min}$ and $T_{max}$ are given in Table 1 respectively.

In the Lactin model (4) a new parametrization is used. Hence, the new
transformed parameters are:

$l=-\lambda\quad,\quad del=\frac{1}{\varDelta}\quad,\quad a={e}^{(\rho-
del)\cdot T_{m}}$

are considered instead of the original parameters $\lambda$, $\varDelta$ and
$T_{m}$, respectively. By their definition, the new parameters are taking
values in $\lambda\in\mathcal{R}$, $0<{del}<1$ and
$a\in\left(0,\rho\cdot\varDelta\right)$. The latter is derived by the
definition of the optimum temperature in the Lactin model (5) and the fact
that $T_{opt}>0$. In addition, should the temperature at the inflection point
$T_{inf}$ be a positive number, then from (6) we can get that
$a\in\left(0,{\rho}^{2}\cdot{\varDelta}^{2}\right)$. The prior distributions
considered are shown in Table 1 respectively for each ecological model.

Table 1: Priors of the parameters of the four ecological models

Models | Parameters | Priors | Models | Parameters | Priors
---|---|---|---|---|---
Bieri | $\alpha$ | U(0,1) | Briere |  †$\widetilde{a}=-\log(\alpha)$ | $\Gamma\left({10}^{-1},{10}^{-2}\right)$
$\beta$ | $\Gamma\left(2\cdot{10}^{-1},{10}^{-1}\right)$ | $T_{min}$ | $\Gamma\left({10}^{-2},{10}^{-2}\right)$
$T_{m_{1}}$ | $\Gamma\left({10}^{-1},{10}^{-2}\right)$ | $T_{max}$ | $\Gamma\left({10}^{-2},{10}^{-3}\right)$
$T_{m_{2}}$ | $\Gamma\left({10}^{-1},{10}^{-2}\right)$ | Analytis |  †$\widetilde{a}=-\log(\alpha)$ | $\Gamma\left({10}^{-1},{10}^{-2}\right)$
Lactin |  †$l=-\lambda$ | $\Gamma\left({10}^{-1},{10}^{-1}\right)$ | $m$ | $\Gamma\left({10}^{-1},{10}^{-1}\right)$
†$del=\frac{1}{\varDelta}$ | U(0,1) | $m$ | $\Gamma\left({10}^{-1},{10}^{-1}\right)$
†$a={e}^{(\rho-del)\cdot T_{m}}$ | $\Gamma\left({10}^{-2},{10}^{-3}\right)$ | $T_{min}$ | $\Gamma\left({10}^{-2},{10}^{-2}\right)$
$\rho$ | U(0,1) | $T_{max}$ | $\Gamma\left({10}^{-2},{10}^{-3}\right)$

  * •

${\dagger}$ transformed parameter used.

##  3 Application and results

We adopt the Bayesian paradigm to inference and model selection. For
concreteness we give the discernible features of the main tools we utilize in
section C in the appendix.

###  3.1 Data application

The reciprocals of the days counted express the observed rates of the insects
and mites from egg to adult stage. The range of the observed rates are at [0,
1] interval. The zero value indicates that no development is observed. This
situation implies the existence of a truncation point which gives explicitly
an upper bound for the upper thermal threshold.

Two datasets are used in this analysis. They concern the study of the two-
spotted mite, Tetranychus urticae (Barber et al., 2003) and the fourteen-
spotted ladybird beetle, Propylea quatuordecimpunctata (Papanikolaou et al.,
2013). The Tetranychus urticae data developmental rates have minimum 0.019 at
$15^{o}C$ and maximum 0.182 at $32.5^{o}C$. The Propylea quatuordecimpunctata
dataset consist of 105 beetles and their developmental rate have minimum 0 at
$35^{o}C$ and maximum 0.111 at $32.5^{o}C$.

###  3.2 Tetranychus urticae

There are 247 mites that have inserted adult stage until the study ended. The
four ecological models are used both assuming the Gaussian and the Inverse
Gamma distributions for the data. The information criteria along with the
estimates of the marginal likelihood for each model are provided in Table 2.

|  | AIC | DIC | LooCV | WAIC | BIC |  $log\left(P_{y}IS\right)^{{\dagger}}$ (se) |  $log\left(P_{y}PP\right)^{{\ddagger}}$ (se) |  $log\left(P_{y}BS\right)^{\S}$ (se)
---|---|---|---|---|---|---|---|---|---
Gaussian model | Bieri | -1663.2 | -1667.7 | -1672.0 | -1672.1 | -1638.6 | 796.1 (28.6) | 777.2 (19.0) | 792.6 (7.1)
Briere | -1592.2 | -1595.8 | -1595.8 | -1595.9 | -1574.7 | 759.6 (39.0) | 728.2 (12.4) | 758.0 (6.9)
Analytis | -1738.3 | -1736.9 | -1744.0 | -1744.0 | -1717.3 | 815.9 (40.5) | 830.2 (11.0) | 821.6 (12.3)
Lactin | -1673.9 | -1675.6 | -1689.1 | -1689.1 | -1638.8 | 797.8 (40.0) | 821.8 (23.2) | 793.8 (56.4)
Inverse Gamma model | Bieri | -1715.6 | -1716.0 | -1717.2 | -1717.2 | -1698.1 | 837.4 (42.6) | 839.7 (9.0) | 832.1 (10.0)
Briere | -1749.4 | -1749.4 | -1746.3 | -1746.3 | -1735.4 | 846.7 (46.1) | 874.0 (12.8) | 848.0 (3.5)
Analytis | -1889.0 | -1887.3 | -1889.0 | -1889.1 | -1867.6 | 904.2 (42.6) | 928.2 (16.3) | 903.9 (26.2)
Lactin | -1911.4 | -1911.0 | -1910.5 | -1910.6 | -1893.9 | 920.9 (30.3) | 956.4 (17.0) | 922.0 (6.5)

  * •

${\dagger}$ $log\left(P_{y}IS\right)$ denotes the logarithm of estimated
marginal likelihood via Importance sampling,

  * •

${\ddagger}$ $log\left(P_{y}PP\right)$ denotes the logarithm of estimated
marginal likelihood via Power posterior,

  * •

$\S$ $log\left(P_{y}BS\right)$ denotes the logarithm of estimated marginal
likelihood via Bridge sampling.

Table 2: Model selection criteria for the eight models applied to the
Tetranychus urticae data

Both Information criteria and marginal likelihood results clearly suggest that
the Inverse Gamma distribution has better fit than the Gaussian distribution
at Tetranychus urticae dataset across all the ecological models. Furthermore
the Lactin model with the Inverse Gamma distribution stands out against all
the other cases. In the Gaussian case, the Analytis model excels.The use of
the Inverse Gamma distribution, on the other hand, not only increases the
Analytis model’s efficiency but also adds flexibility to the Lactin model. The
Briere model has the poorest criteria values (Table 2). Although the
information criteria have fairly indicated the Briere model as the most
suitable for the data, when we concentrate on the upper thermal threshold
$T_{max}$ estimate, its performance is poor compared to the other ecological
models shown in Fig. 5. Nevertheless, even though there is a clear picture
concerning information criteria values between the ecological models, there is
some variation between marginal likelihood estimates within ecological models
in Table 2.

Posterior means, $95\%$ credible limits and the effective sample size (neff)
of the thermal thresholds and the deviance of the four ecological models are
summarised in Tables 3.2 and A.1 using the Inverse Gamma and the the Gaussian
distribution, respectively. In addition, the HMC, ADVI meanfield and ADVI
fullrank methods appear alternately in each column for each parameter of
interest.

The $T_{min}$ credible limits estimates between the four ecological models do
not overlap in Tables 3.2 and A.1. Bieri model has greater limits and Lactin
model gives negative value estimates in the Inverse Gamma case Table 3.2. The
credible limits for $T_{opt}$ overlap between Bieri and Analytis whereas the
estimates are lower for Lactin model and greater for Briere model. The
$T_{max}$ credible limits overlap for Bieri and Lactin in Table A.1, has
higher values for Briere and lower for Analytis model.

The credible limits using ADVI meanfield and fullrank methods overlap with the
HMC credible limits in most cases in Tables 3.2 and A.1. Also the ADVI
fullrank method seems to be closer to the HMC estimates as in Briere and
Lactin models in Table 3.2. However in general it has worst fit than the
corresponding fitted model using HMC and also gives wider $95\%$ credible
intervals.

Bayesian model averaging provides alternative estimates for the parameters of
interest, combining the predictive efficiency of all four ecological models.
The derived weights are shown in Table A.1 whereas the BMA estimates and their
$95\%$ credible intervals are given on Table A.1 and are divided into the data
case of the Gaussian distribution and the data case of the Inverse Gamma
distribution. The predictive bias of the Analytis and Lactin models appear to
affect the model averaging estimates of $T_{min}$, $T_{max}$ and the deviance
in the Gaussian and the Inverse Gamma case respectively. Also different
weights based on (20) give almost identical $95\%$ credible limits. However
using ELBO based BMA weights do not give robust estimates of the parameters of
interest which is not unexpected since the ELBO is a lower bound estimate of
the marginal likelihood. Furthermore, the time elapsed until the completion of
the algorithm for ADVI methods is up to 52 seconds, while for the HMC method
is at least 476 seconds for Tetranychus dataset as shown in Table A.1. The
average difference between the working times of HMC and each ADVI method is
around 3500 seconds (58 minutes) whereas between the meanfiled and fullrank is
around 18 seconds. These findings illuminate the VBI fullrank method’s time-
efficacy.

Figure 5: Predicted posteriors versus Tetranychus urticae data using Inverse
Gamma distribution

Posterior summaries for the four models using the Inverse Gamma distribution
for the Tetranychus urticae data. In each column we report the HMC, the ADVI-
meanfield and ADVI-fullrank estimates respectively. $T_{min}$ neff‡ $T_{opt}$
neff‡ Mean & $95\%$ Cr. I. Bieri 9.4 (9.3, 9.6) 9.3 (9.1, 9.4) 9.3 (9.0, 9.5)
19741 33.7 (33.2, 34.1) 145.5 (134.8, 149.8) 203.2 (114.5, 326.6) 8888 Briere
6.6 (6, 7) 6.6 (6.4, 6.7) 6.5 (6.1, 7) 7933 36.7 (35.2, 38.5) 36.6 (36.1, 37)
36.6 (35, 38.3) 6978 Analytis 4.2 (4, 4.6) 7.6 (7.5, 7.8) 6.0 (4.4, 9.4) 9977
33.6 (33.3, 34) 99.9 (96.6, 103) 79.9 (16.8, 239.3) 7876 Lactin -18.7 (-18.7,
-18.7) 7.8 (7.6, 8.1) -18.7 (-18.8, -18.7) 18393 32.0 (31.8, 32.2) 32.1 (31.8,
32.3) 33.9 (33.8, 34.1) 12228 $T_{max}$ neff‡ dev† neff‡ Mean & $95\%$ Cr. I.
Bieri 35.8 (35.4, 36.5) 150.8 (149.2, 152.5) 209.0 (114.5, 326.6) 7738 -1720.9
(-1724.9, -1713) -1652.9 (-1658.6, -1637.6) -1652.4 (-1658.5, -16) 10228
Briere 45.0 (43.1, 47.3) 44.8 (44.2, 45.4) 44.9 (42.8, 47) 6935 -1753.3
(-1757, -1745.9) -1752.9 (-1757.1, -1741.4) -1751.9 (-1756.7, -1740.4) 8638
Analytis 35.0 (35, 35.1) 99.9 (96.6, 103.1) 80.9 (17.5, 240.6) 9038 -1894.0
(-1899.2, -1885) -1670.9 (-1691.8, -1641.4) -1495.1 (-1687.4, -149.5) 12370
Lactin 38.4 (38.1, 38.9) 42.6 (42.4, 42.8) 38.6 (38.2, 39) 8284 -1916.3
(-1920.6, -1908.3) -1705.8 (-1747.8, -1650.5) -1890.0 (-1918.5, -1791.8) 9938
• ${\dagger}$ deviance of the model given the data. • ${\ddagger}$ effective
sample size.

In Fig. 5 the predicted posteriors versus Tetranychus urticae data using
Inverse Gamma distribution are shown for the ecological models.

The adaptivity in data across predictor values is clear in all the ecological
models, where the $95\%$ credible limits are adjusted to data variance in each
temperature level. On the other hand, in the Gaussian case (Fig. 7) the
variance of the posterior remains constant across predictor values.

###  3.3 Propylea quatuordecimpunctata dataset

There are 17 out of 105 insects that have not altered their egg status until
the study ended. These cases are observed at ${35^{\circ}C}$ and are indicated
by zeros in the response variable y. The Inverse Gamma distribution is not
defined for zero response values. On the other hand, the Bieri and Lactin
models can also generate negative values. Hence, in order to account for the
presence of zeros, we use a Zero Inflated Inverse Gamma model and choose
parameter initial values so as the scale of the Inverse Gamma to remain
positive. Both the Gaussian and the Zero Inflated Inverse Gamma distributions
are used for the data along with the four ecological models. The various model
selection tools are summarised in Table A1.

Both Information criteria and marginal likelihood results clearly suggest that
the Zero Inflated Inverse Gamma distribution has better fit than the Gaussian
distribution at Propylea quatordicempuncata only at the Analytis model.
Furthermore Bieri and Analytis models stand out in the Gaussian and the
Inverse Gamma case respectively, whilst the Briere model has the poorest
criteria values (Table A1). Nevertheless, even though there is a clear picture
of the goodness of fit between the ecological models, there is some variation
between marginal likelihood estimates within ecological models in Table A1.

In the Gaussian model, the Bieri and Lactin models stand out according to the
information criteria and marginal likelihood estimates, the Analytis model
follows, whilst the Briere model has the lower criteria values. On the
contrary, in the Zero Inflated Inverse Gaussian model, the Analytis model is a
better choice according to the marginal likelihood and the information
criteria values, while Bieri, Lactin and Briere models follow respectively.
Lactin and Bieri can be interpreted by their ability to track an exponential
data-value decrease. Additionally, the Bieri model manages to capture the
linear increase in the Propylea quatuordecimpunctata dataset. The Briere
model, on the other hand, lacks performance due to its unique structure, which
requires that the decline of the response variable as the temperature rises be
of the form $\sqrt{T_{max}-T}$. The Analytis model has similar multiplicative
structure to the Briere model, but it is more complex model since the
exponents of its model are unknown variables that make it adaptive and liable
to data change especially in the Inverse Gamma case. Nevertheless, the
Analytis model needs some tuning when defining its thermal parameters space.
This is necessary to avoid allowing $T_{min}$ values close to zero, which
would result in parameter underestimation and poor fit in the Analytis model.
Especially, in the case of the Propylea dataset, we a-prior assumed that
$T_{min}$ is greater than ${4^{\circ}C}$. Posterior means, $95\%$ credible
limits and the effective sample size of the thermal thresholds and the
deviance of the four ecological models are summarised in Tables 3.3 and A.1
using the Zero inflated Inverse Gamma and the Gaussian distribution
respectively. Specifically, the HMC, ADVI meanfield and ADVI fullrank methods
appear alternatively in each column of Tables (3.3 and A.1) for each parameter
of interest. Briere model has greater 95% credible limits for $T_{min}$ while
Lactin model gives negative value estimates. The results for $T_{opt}$ overlap
for Bieri and Analytis in the Zero Inflated Inverse Gamma case whereas the
estimates are lower for Lactin model. The $T_{max}$ credible limits overlap
between Bieri and the other ecological models both in the Gaussian and the
Zero Inflated Inverse Gamma case. The zero-rate values can be naturally
modelled in the case of the Gaussian distribution. Though only the Bieri and
Analytis models give wider $95\%$ credible intervals. For the rest of the
models, the credible limits difference is lower than three degrees
${3^{\circ}C}$. In addition, in Zero Inflated Inverse Gamma case, the Analytis
model has higher information criteria values, while the Bieri and Lactin
models follow as shown in Table A1, respectively. All four ecologocal models
incorporate terms into their structure in the form of products which directly
affect the scale of the Inverse Gamma distribution. Thus, not only the mean
but also the variance of the statistical model change along and adapt the
flunctuations for each temperature level as shown in Fig. 9. This adaptivity
is evident across all the four ecological models in use.

Figure 6: Boxplots of the posterior probability of non-zero development in
Propylea Coccinellidae data using the Zero Inflated Inverse Gamma model

The $T_{min}$ outcomes for 95% credible limits do not overlap, while Briere
model has larger values. The Lactin model is not suitable to provide
reasonable estimates as it allows for negative values (Table 3.3). The results
for $T_{opt}$ are sorted in ascending order from Bieri, Lactin, Analytis, and
Briere. We observe that the credible limits for the $T_{max}$ overlap between
models and do not exceed the maximum observed temperature $35^{o}C$.
Concerning the Zero Inflated Inverse Gamma case, the posterior predictive
probabilities of non-zero entries, which are the probabilities of anthropods’
development, are shown in Figure 6. We observe that the probabilities tend to
unity near $T_{opt}$ estimates whereas become negligible out of the thermal
threshold limits except for the Lactin model case which does not give robust
estimates for $T_{min}$. In addition, the 95% credible limits estimates using
ADVI meanfield and fullrank methods overlap with the HMC estimates with
exceptions at the Lactin model in which the fullarank method seems to be
closer to the HMC estimates (Tables A.1 and 3.3). In many cases, the ADVI
fullrank method agrees with HMC. It seems to give more robust estimates for
$T_{max}$ mean in Gaussian case and Bieri model. However it has worst fit than
HMC and also gives wider $95\%$ credible intervals.

Posterior summaries for the four models using the Zero Inflated Inverse Gamma
distribution for the Propylea Coccinellidae data. In each column we report the
HMC, the ADVI-Mean field and ADVI-Full rank estimates respectively. $T_{min}$
neff‡ $T_{opt}$ neff‡ Mean & 95% Cr. I. Bieri 9.9 (9.4, 10.3) 9.4 (9.3, 9.6)
9.2 (8.7, 9.7) 9138 32.7 (32.2, 33.8) 73.0 32.6, 146.6 109.4 (28.7, 320.1)
16552 Briere 11.1 (10.3, 11.9) 11.2 (10.8, 11.5) 11.2 (10.4, 12.0) 14192 29.3
(29.2, 29.5) 29.3 (29.2, 29.4) 29.3 (29.2, 29.5) 17300 Analytis 5.0 (4.0, 7.0)
5.1 (4.1, 8.7) 4.4 (4.2, 4.9) 352 33.5 (32.3, 34.9) 33.0 (29.71, 37.0) 32.2
(30.5, 33.9) 1240 Lactin -133.51 (-157.5, -116.8) 11.6 (1.2, 59.2) -12.5
(-409.8, 8.5) 11511 30.9 (30.8, 31.1) 39.7 (1.1, 211.2) 25.9 (5.5, 77.7) 9847
$T_{max}$ neff‡ dev† neff‡ Mean & 95% Cr. I. Bieri 34.4 (33.7, 34.9) 18.1
(9.3, 52.6) 16.0 (8.8, 52.5) 12972 -722.3 (-725.7, -714.8) -318.5 (-360.4,
-318.8) -323.8 (-572.1, -309.3) 9962 Briere 35.0 (34.8, 35.0) 34.9 (34.8,
35.0) 34.8 (34.8, 35.0) 23972 -563.8 (-568.0, -555.9) -563.1 (-567.7, -552.0)
-561.9 (-567.4, -549.7) 14077 Analytis 33.6 (32.5, 34.9) 33.6 (30.3, 37.6)
32.5 (30.7, 34.4) 1409 -736.2 (-742.1, -729.1) -411.6 (-482.3, -216.9) -684.5
(-742.5, -618.8) 1117 Lactin 34.9 (34.7, 35.0) 25.3 (1.1, 76.3) -45.1 (-935.4,
52.5) 16361 -703.6 (-708.1, -695.0) 247.8 (0.0, 294.8) 247.8 (0.0, 294.8) 9179
• ${\dagger}$ deviance of the model given the data, • ${\ddagger}$ effective
sample size.

Furthermore, the time elapsed until the completion of the algorithm for ADVI
methods is up to 48 seconds, while for the HMC method it is at least 321
seconds for Propylea dataset as shown in Table A.1. The average difference
between the working times of HMS and each ADVI method is around 6300 seconds
(105 minutes) whereas between the meanfield and fullrank is around 12 seconds.
These findings illuminate both the VBI fullrank method’s computational time-
efficacy.

Bayesian model averaging results are shown in Tables A.1, A.1. This time the
predictive bias that affect the model averaging estimates are from the Bieri
and Analytis models fot he Gaussian and the Zero Inflated Inverse Gamma
distribution respectively.

##  4 Conclusions

Conclusions of the current work can be summarised into five groups.

###  4.1 Comparison of Computational Methods

Although ADVI methods agree with HMC in many of cases and are distinguished in
time-efficacy, robust estimates do not appear to be available for all
parameters of interest, with the ADVI fullrank method providing better
estimates than the ADVI meanfield mehod. The two datasets application show
that they are sensitive to the initialization of the parameters and the root-
finding algorithm used, in particular for more complicated models such as
Analytis and Lactin. For example, if the root initial suggestion is positive,
the root generated is also positive, even if the HMC results are negative,
like in Tables(3.2, 3.3, A.1). However, in the same tables, the VBI fullrank
method’s credible limits overlap with those of the HMC method more frequently
than the VBI meanfield method, indicating that the ADVI fullrank method has
little better computational efficiency.

###  4.2 Distribution of the data

Inverse Gamma distribution is a promising alternative to the classical
Gaussian distribution in the analysis of the rate of development of anthropods
in the absence of zero data and also provides adaptive temperature-level
variance predictions across all ecological models. It also reinforces the
flexibility of the Lactin model. In addition Inverse Gamma distribution
naturally models the developmental rates that are defined as the reciprocal of
positive real numbers.

###  4.3 Zero-inflation performance

Zero Inflated Inverse Gamma distribution can be used in the presence of non-
development (zero data) after defining the probability of generating zeros. In
this way the developmental rates are naturally modeled and also the
predictions are variance adaptive to the temperature levels. Finally, it
performs better than the classical Gaussian model in the Analytis case as far
as information criteria scores and marginal likelihood estimates are
concerned.

###  4.4 Model comparison

For the model comparison, both the information criteria and the marginal
probability values agree on the prioritization of performance across
ecological models at both data-sets. Nevertheless, there is some variation
between marginal likelihood estimates within ecological models and their is no
clear pattern across the marginal likelihood approaches.

There is no standard ecological model that has better fit measured by
information criteria and the estimated marginal likelihood than the others in
all the situations. However, even in the presence of zero response values, the
Bieri ecological model appears to have robust results using the Gaussian
distribution and stands out when there is a linear increase of the
developmental rate like in the Propylea quatuordecimpunctata data-set case. In
addition, it provides interpretable estimates of interest parameters with
credible intervals, often overlapping with that of the Analytis or the Lactin
model, using weakly informative priors. By contrast, the Bieri model can not
follow the non-linear increase in the rate of development of the Tetranychus
data, as do the more sophisticated Analytis or / and Lactin models.
Lactin model has better fit in the case of exponentially fluctuations of the
developmental rate, but is inadequate for $T_{min}$ estimation. Analytis model
needs tuning for $T_{min}$ sample space so as not to underestimate it and
provides robust estimates for the other thermal parameters involved in the
model.

Briere and Analytis models by definition includes the true threshold
parameters $T_{min}$ and $T_{max}$ along with a functional truncation between
these threshold parameters, which do not affect their performance. Moreover,
the Briere model appears to have the worst marginal likelihood values, as well
as the worst predictive performance for both Gaussian and Inverse Gamma
distributions, despite having the least parameters and being widely used. It
also seems to generate estimates with small 95% credible intervals and higher
$T_{max}$ values. Although statistical inference based on the information
criteria have reasonable values for Briere model, when we focus on a
particular parameter of interest such as the $T_{max}$, its performance may be
poor compared to the other ecological models as shown in Fig. 5. In
conclusion, the Briere model can be used to obtain quick estimates of thermal
parameters involved in arthropod evolution without the need for any pre-
adjustment (parameters tuning or / and transformation).

The Analytis model needs some initial boost by either tuning the thermal
parameters or using informative priors in order to achieve good predictions.
In particular, a lower bound of $4^{o}C$ for $T_{min}$ should be imposed in
order not to generate negligible or negative values. An upper bound for
$T_{max}$ should also be imposed if there are zero responses at a higher
temperature. It is only after these parameters have been calibrated that the
Analytis model stands out in its performance by comparing the other models and
generates wide 95% credible intervals for parameters of interest.
Additionally, due to its functional form, Analytis model can fit any
polynomial shaped increase and decrease of the developmental rate. In
comparison to the simple Gaussian case, Analytis explains the zeros generated
scheme better. However, it provides smaller estimates than the other three
models, for the $T_{min}$ parameter. In summary, the Analytis model performs
well to predict the upper thermal parameters only after adjusting the
$T_{min}$ and $T_{max}$ sample space interval limits to prior information from
the data.

In addition, the Lactin model does not include an analytical lower thermal
threshold in its form. The lower thermal threshold estimates are implicitly
derived as the lower root of the Lactin model applying the appropriate model
parameters values. Thus, the Lactin model in some cases generates negative
estimates of $T_{min}$ that can not be interpreted in ecology studies. Also
parameters $\lambda$ and $\varDelta$ fluctuations trigger statistical model
discontinuity situations that need to be addressed taking into account the
data-set structure. Moreover, it has smaller effective sample size especially
in the case of Gaussian distribution. On the contrary, it generates a wider
95% credible interval for $T_{max}$ compared to other ecological models and
also performs better in the case of zero responses at high temperatures. In
summary, the Lactin model is applicable when research focus on robust
estimates of the optimum and the upper thermal threshold parameters without
the need for prior tuning.

###  4.5 BMA performance

BMA technique is used to address the difficulty that the true causal structure
in ecological data is often not known and often several models are used to
describe the development of anthropods. The mean square error in model-
averaged predictions depends on each model’s predictive bias and variance, the
covariance in predictions between models, and the uncertainty about model
weights (Dormann et al., 2018). In this work, we adopt weights in (20). The
BMA weights summarised in Tables(A.1, A.1) clearly support the Bieri and
Analytis models for the Propylea data and the Analytis and Lactin models for
the Tetranychus data. As a result, the BMA estimates are close to the
estimates of these models respectively for the two data-sets.

##  5 Discussion

There are several issues that we deal with in terms of arthropod developmental
rates, and there are some concerns for future research. To begin with,
comparing non-linear non-nested models with varying numbers of parameters and
truncated mean structures, as well as excessive-zeros in data, is not a
trivial task.

###  5.1 Comparison of Computational Methods

We use the Bayesian paradigm, as well as some contemporary computational
approaches such as the HMC, ADVI-meanfield, and ADVI-fullrank, to address not
only irregular and truncated mean structures, but also the uncertainty of zero
generation in the data. Although ADVI techniques are gaining popularity in the
scientific community due to their fast and computationally inexpensive
approximations to the posterior distributions, they do not provide robust
estimates of all the parameters of the models we study. The HMC method, on the
other hand, gives robust estimates even under these specific model and data
structure conditions.

###  5.2 Distribution of the data

Furthermore, for the data generation scheme, we suggest the Gaussian and the
Inverse Gamma distributions. The Gaussian option provides sensible estimates
that can be used even though the data has a lot of zeros. Inverse Gamma, on
the other hand, not only naturally models developmental rates, which are
characterized as the reciprocal of positive real values, but also provides
variance adaptivity across temperature fluctuations. Allso for each ecological
model we define the Zero Inflated Inverse Gamma density so as to model data
with an excessive number of zeros. When comparing models involving Gaussian
and Inverse Gamma or Zero Inflated gamma distributions, we find that the
second performs better in non-zero data cases, while the first performs better
in all but the Analytis model.

###  5.3 Model comparison

In addition, we address the model comparison challenge by employing the
Information criteria not only to assess model goodness of fit to data while
accounting for model complexity, but also to assess ability of models to make
predictions on new data using the leave one out cross validation technique.
Additionally, we use marginal likelihood approximations of the various models
to determine which one is best supported by the data. Finally, we plot the
predicted posteriors alongside the observed data points to visualize the
prediction ability of the suggested models.

###  5.4 BMA performance

The predictive bias of the weighted models, as well as the uncertainty about
the weights, affect BMA results using Information criteria as weights
according to (20) weights as outlined in section C.2.3 in the appendix. As a
result, the BMA approach does not provide estimates that differ from the best-
performing models.

###  5.5 Future research

Among things for future research is to select consistently the most robust
candidate between models given sufficiently many data samples, in a
sensitivity analysis perspective. Moreover, the ADVI methods, can be extended
so as to capture more sophisticated mean structures, like the ones we present
in the current work. In addition, probability density that generates zeros in
the Zero Inflated Inverse Gamma distribution can be modeled in more complex
ways, such as using hyperparameters and hierarchical effects across
temperature levels. Finally, R-packages that include the suggested models and
perform the analysis presented in this paper are to be created.

## References

  * Aghdam et al. (2011) Aghdam, H. R., Fathipour, Y. and Kontodimas, D. (2011) Evaluation of non-linear models to describe development and fertility of codling moth at constant temperatures.  Entomologia Hellenica, 20, 3–16.
  * Akaike (1974) Akaike, H. (1974) A new look at the statistical model identification.  IEEE Transactions on Automatic Control, 19, 716–723.
  * Analytis (1981) Analytis, S. (1981) Relationship between temperature and development times in phytopathogenic fungus and in plant pests: a mathematical model.  Agric. Res. Athens, 5, 133–159.
  * Barber et al. (2003) Barber, A., Campbell, C., Crane, H., Lilley, R. and Tregidga, E. (2003) Biocontrol of two-spotted spider mite tetranychus urticae on dwarf hops by the phytoseiid mites phytoseiulus persimilis and neoseiulus californicus.  Biocontrol Science and Technology, 13, 275–284.
  * Betancourt (2017) Betancourt, M. (2017) A conceptual introduction to hamiltonian monte carlo.  arXiv preprint arXiv:1701.02434.
  * Bieri et al. (1983) Bieri, M., Baumgartner, J., Bianchi, G., Delucchi, V., Arx, R. v. et al. (1983) Development and fecundity of pea aphid (acyrthosiphon pisum harris) as affected by constant temperatures and by pea varieties.  Mitteilungen der Schweizerischen Entomologischen Gesellschaft, 56, 163–171.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A. and McAuliffe, J. D. (2017) Variational inference: A review for statisticians.  Journal of the American Statistical Association, 112, 859–877.
  * Bradshaw et al. (2016) Bradshaw, C. J., Leroy, B., Bellard, C., Roiz, D., Albert, C., Fournier, A., Barbet-Massin, M., Salles, J.-M., Simard, F. and Courchamp, F. (2016) Massive yet grossly underestimated global costs of invasive insects.  Nature communications, 7, 12986.
  * Briere et al. (1998) Briere, J., Pracros, P. and Stoeckel, J. (1998) Modeling development rate for predicting lobesia botrana (den. & schiff.) population dynamics.  IOBC WPRS Bulletin, 21, 51–52.
  * Briere and Pracros (1998) Briere, J.-F. and Pracros, P. (1998) Comparison of temperature-dependent growth models with the development of lobesia botrana (lepidoptera: Tortricidae).  Environmental Entomology, 27, 94–101.
  * Briere et al. (1999) Briere, J.-F., Pracros, P., Le Roux, A.-Y. and Pierre, J.-S. (1999) A novel rate model of temperature-dependent development for arthropods.  Environmental Entomology, 28, 22–29.
  * Broufas et al. (2007) Broufas, G., Pappas, M. and Koveos, D. (2007) Development, survival, and reproduction of the predatory mite kampimodromus aberrans (acari: Phytoseiidae) at different constant temperatures.  Environmental entomology, 36, 657–665.
  * Broufas and Koveos (2001) Broufas, G. D. and Koveos, D. S. (2001) Development, survival and reproduction of euseius finlandicus (acari: Phytoseiidae) at different constant temperatures.  Experimental & applied acarology, 25, 441–460.
  * Brown et al. (2004) Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. and West, G. B. (2004) Toward a metabolic theory of ecology.  Ecology, 85, 1771–1789.
  * Buckland et al. (1997) Buckland, S. T., Burnham, K. P. and Augustin, N. H. (1997) Model selection: an integral part of inference.  Biometrics, 603–618.
  * Campbell et al. (1974) Campbell, A., Frazer, B., Gilbert, N., Gutierrez, A. and Mackauer, M. (1974) Temperature requirements of some aphids and their parasites.  Journal of Applied Ecology, 11, 431–438.
  * Carlin and Chib (1995) Carlin, B. P. and Chib, S. (1995) Bayesian model choice via markov chain monte carlo methods.  Journal of the Royal Statistical Society: Series B (Methodological), 57, 473–484.
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P. and Riddell, A. (2017) Stan: A probabilistic programming language.  Journal of statistical software, 76, 1–32.
  * Damos and Savopoulou-Soultani (2012) Damos, P. and Savopoulou-Soultani, M. (2012) Temperature-driven models for insect development and vital thermal requirements.  Psyche: A Journal of Entomology, 2012, 13.
  * Dellaportas et al. (2000) Dellaportas, P., Forster, J. J. and Ntzoufras, I. (2000) Bayesian variable selection using the gibbs sampler.  Biostatistics-Basel, 5, 273–286.
  * Dormann et al. (2018) Dormann, C. F., Calabrese, J. M., Guillera-Arroita, G., Matechou, E., Bahn, V., Bartoń, K., Beale, C. M., Ciuti, S., Elith, J., Gerstner, K. et al. (2018) Model averaging in ecology: a review of bayesian, information-theoretic, and tactical approaches for predictive inference.  Ecological Monographs, 88, 485–504.
  * Friel and Pettitt (2008) Friel, N. and Pettitt, A. N. (2008) Marginal likelihood estimation via power posteriors.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70, 589–607.
  * Frühwirth-Schnatter (2004) Frühwirth-Schnatter, S. (2004) Estimating marginal likelihoods for mixture and markov switching models using bridge sampling techniques.  The Econometrics Journal, 7, 143–167.
  * Gelman et al. (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A. and Rubin, D. B. (2013) Bayesian data analysis.  CRC press.
  * Gelman and Meng (1998) Gelman, A. and Meng, X.-L. (1998) Simulating normalizing constants: From importance sampling to bridge sampling to path sampling.  Statistical science, 13, 163–185.
  * George and McCulloch (1997) George, E. I. and McCulloch, R. E. (1997) Approaches for bayesian variable selection.  Statistica sinica, 7, 339–373.
  * Green (1995) Green, P. J. (1995) Reversible jump markov chain monte carlo computation and bayesian model determination.  Biometrika, 82, 711–732.
  * Gronau et al. (2020) Gronau, Q. F., Singmann, H. and Wagenmakers, E.-J. (2020) bridgesampling: An R package for estimating normalizing constants.  Journal of Statistical Software, 92, 1–29.
  * Hare et al. (2011) Hare, M. P., Nunney, L., Schwartz, M. K., Ruzzante, D. E., Burford, M., Waples, R. S., Ruegg, K. and Palstra, F. (2011) Understanding and estimating effective population size for practical application in marine species management.  Conservation Biology, 25, 438–449.
  * van der Have (2008) van der Have, T. M. (2008) Slaves to the Eyring equation?: Temperature dependence of life-history characters in developing ectotherms.  Ph.D. thesis, Wageningen University.
  * Huey and Berrigan (2001) Huey, R. B. and Berrigan, D. (2001) Temperature, demography, and ectotherm fitness.  The American Naturalist, 158, 204–210.
  * Jarošík et al. (2002) Jarošík, V., Honěk, A. and Dixon, A. F. (2002) Developmental rate isomorphy in insects and mites.  The American Naturalist, 160, 497–510.
  * Kass and Raftery (1995) Kass, R. E. and Raftery, A. E. (1995) Bayes factors.  Journal of the american statistical association, 90, 773–795.
  * Kontodimas et al. (2004) Kontodimas, D. C., Eliopoulos, P. A., Stathas, G. J. and Economou, L. P. (2004) Comparative temperature-dependent development of nephus includens (kirsch) and nephus bisignatus (boheman)(coleoptera: Coccinellidae) preying on planococcus citri (risso)(homoptera: Pseudococcidae): evaluation of a linear and various nonlinear models using specific criteria.  Environmental Entomology, 33, 1–11.
  * Kuo and Mallick (1998) Kuo, L. and Mallick, B. (1998) Variable selection for regression models.  Sankhyā: The Indian Journal of Statistics, Series B, 60, 65–81.
  * Lactin et al. (1995) Lactin, D. J., Holliday, N., Johnson, D. and Craigen, R. (1995) Improved rate model of temperature-dependent development by arthropods.  Environmental entomology, 24, 68–75.
  * Meng and Wong (1996) Meng, X.-L. and Wong, W. H. (1996) Simulating ratios of normalizing constants via a simple identity: a theoretical exploration.  Statistica Sinica, 6, 831–860.
  * Neal (2010) Neal, R. M. (2010) MCMC using Hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 54, 113–162.
  * Nedvěd (2009) Nedvěd, O. (2009) Temperature, effects on development and growth.  In Encyclopedia of Insects (Second Edition), 990–993. Academic Press, second edition edn.  URL: https://www.sciencedirect.com/science/article/pii/B9780123741448002617.
  * Neven (2000) Neven, L. G. (2000) Physiological responses of insects to heat.  Postharvest Biology and Technology, 21, 103–111.
  * Norris and Kunz (2012) Norris, A. L. and Kunz, T. H. (2012) Effects of solar radiation on animal thermoregulation.  Solar radiation, 195–220.
  * Overstall and Forster (2010) Overstall, A. M. and Forster, J. J. (2010) Default bayesian model determination methods for generalised linear mixed models.  Computational Statistics & Data Analysis, 54, 3269–3288.
  * Papanikolaou et al. (2019) Papanikolaou, N. E., Kavallieratos, N. G., Kondakis, M., Boukouvala, M. C., Nika, E. P. and Demiris, N. (2019) Elucidating fitness components of the invasive dermestid beetle trogoderma granarium combining deterministic and stochastic demography.  PloS one, 14, e0212182.  URL: https://doi.org/10.1371/journal.pone.0212182.
  * Papanikolaou et al. (2013) Papanikolaou, N. E., Milonas, P. G., Kontodimas, D. C., Demiris, N. and Matsinos, Y. G. (2013) Temperature-dependent development, survival, longevity, and fecundity of propylea quatuordecimpunctata (coleoptera: Coccinellidae).  Annals of the Entomological Society of America, 106, 228–234.
  * Papanikolaou et al. (2014) Papanikolaou, N. E., Milonas, P. G., Kontodimas, D. C., Demiris, N. and Matsinos, Y. G. (2014) Life table analysis of propylea quatuordecimpunctata (coleoptera: Coccinellidae) at constant temperatures.  Annals of the Entomological Society of America, 107, 158–162.
  * Perrakis et al. (2014) Perrakis, K., Ntzoufras, I. and Tsionas, E. G. (2014) On the use of marginal posteriors in marginal likelihood estimation via importance sampling.  Computational Statistics & Data Analysis, 77, 54–69.
  * R Core Team (2020) R Core Team (2020) R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.  URL: https://www.R-project.org/.
  * Schwarz et al. (1978) Schwarz, G. et al. (1978) Estimating the dimension of a model.  The annals of statistics, 6, 461–464.
  * Spiegelhalter et al. (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and Van Der Linde, A. (2002) Bayesian measures of model complexity and fit.  Journal of the royal statistical society: Series b (statistical methodology), 64, 583–639.
  * Wagner et al. (1991) Wagner, T., Olson, R. and Willers, J. (1991) Modeling arthropod development time.  Journal of Agricultural Entomology, 8, 251–270.
  * Watanabe (2010) Watanabe, S. (2010) Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory.  Journal of Machine Learning Research, 11, 3571–3594.

##  Appendix A Tables with summaries and posterior estimates of the
statistical methods used

###  A.1 Propylea quatuordecimpunctata dataset

Algorithmic working time in seconds, for two datasets. Tetranychus urticae
Gaussian Bieri Briere Analytis Lactin HMC 2011 476 9852 5672 ADVI-meanfield 8
1 52 1 ADVI-fullrank 45 6 38 44 Inverse Gamma HMC 1375 761 3231 4917 ADVI-
meanfield 2 4 15 27 ADVI-fullrank 2 34 14 44 Propylea quatuordecimpunctata
Gaussian Bieri Briere Analytis Lactin HMC 4335 20218 4913 1426 ADVI-meanfield
13 7 3 16 ADVI-fullrank 48 11 12 16 Zero Inflated Inverse Gamma HMC 1078 321
15291 3073 ADVI-meanfield 41 5 9 10 ADVI-fullrank 52 5 39 18

BMA weights for the Tetranychus urticae data Gaussian distribution aic_w dic_w
loocv_w waic_w bic_w elbo_mf elbo_fr Bieri 4.9E-17 9.4E-16 2.3E-16 2.4E-16
8.1E-18 3.1E-53 7.3E-45 Briere 1.9E-32 2.3E-31 6.6E-33 6.9E-33 1.1E-31 2.9E-27
1.5E-18 Analytis 1 1 1 1 1 0.001 1 Lactin 1.0E-14 4.9E-14 1.2E-12 1.2E-12
9.0E-18 0.999 3.5E-39 Inverse Gamma distribution Bieri 3.0E-43 4.5E-43 1.1E-42
1.0E-42 3.0E-43 1.6E-20 1.0E-200 Briere 6.6E-36 8.1E-36 2.2E-36 2.1E-36
3.8E-35 8.9E-19 1.0E-171 Analytis 1.4E-5 7.1E-6 2.1E-5 2.1E-5 1.9E-6 1 1
Lactin 0.999 0.999 0.999 0.999 1 1.0E-36 2.7E-58

BMA weights for the Propylea quatuordecimpunctata data Gaussian distribution
aic_w dic_w loocv_w waic_w bic_w elbo_mf elbo_fr Bieri 0.999 0.968 0.991 0.991
1 1 1 Briere 3.3E-34 7.3E-30 5.5E-35 5.5E-35 4.8E-33 5E-203 3E-218 Analytis
1.1E-31 1.5E-32 3.8E-32 3.8E-32 4.0E-31 9.0E-269 1.0E-210 Lactin 5.9E-04 0.032
0.009 0.009 1.1E-05 5.0E-159 3.0E-162 Zero Inflated Inverse Gamma distribution
Bieri 9.3E-5 0.001 0.001 0.001 3.5E-4 5.1E-53 3.4E-72 Briere 1.7E-38 3.9E-38
7.3E-38 7.2E-38 2.4E-37 1 1.1E-23 Analytis 1 0.999 0.999 0.999 1 3.3E-35 1
Lactin 1.7E-8 2.1E-8 5.8E-8 5.8E-8 6.4E-8 3.0E-114 2.0E-140

The mean and $95\%$ Cr.I. limits of the BMA estimates of parameters of
interest calculated with the use of information criteria score weights
assuming the Gauussian and Inverse Gamma distributions respectively for
Tetranychus urticae data. Gaussian Model aic_w dic_w loocv_w waic_w bic_w
elbo_mf elbo_fr $T_{min}$ Mean 4.4 4.4 4.4 4.4 4.4 4.4 4.4 $95\%$ Cr.I. (4.0,
5.3) (4.0, 5.3) (4.0, 5.3) (4.0, 5.3) (4.0, 5.3) (4.0, 5.3( (4.0, 5.3)
$T_{opt}$ Mean 32.9 32.9 32.9 32.9 32.9 32.6 32.9 $95\%$ Cr.I. (32.7, 33.3)
(32.7, 33.3) (32.7, 33.3) (32.7, 33.3) (32.7, 33.3) (32.4, 32.8) (32.7, 33.3)
$T_{max}$ Mean 35.3 35.3 35.3 35.3 35.3 36.9 35.3 $95\%$ Cr.I. (35.1, 35.5)
(35.1, 35.5) (35.1, 35.5) (35.1, 35.5) (35.1, 35.5) (36.5, 37.3) (35.1, 35.5)
dev† Mean -1750.3 -1750.3 -1750.3 -1750.3 -1750.3 -1693.9 -1750.3 $95\%$ Cr.I.
(-1758.0,-1738.3) (-1758.0,-1738.3) (-1758.0,-1738.3) (-1758.0,-1738.3)
(-1758.0,-1738.3) (-1703.1,-1679.7) (-1758.0,-1738.3) Inverse Gamma Model
$T_{min}$ Mean 4.2 4.2 4.2 4.2 4.2 4.2 4.2 $95\%$ Cr.I. (4.0, 4.6) (4.0, 4.6)
(4.0, 4.6) (4.0, 4.6) (4.0, 4.6) (4.0, 4.6) (4.0, 4.6) $T_{opt}$ Mean T_opt
32.0 32.0 32.0 32.0 32.0 33.6 33.6 $95\%$ Cr.I. (31.8, 32.2) (31.8, 32.2)
(31.8, 32.2) (31.8, 32.2) (31.8, 32.2) (33.3, 33.9) (33.3, 33.9) $T_{max}$
Mean 38.4 38.4 38.4 38.4 38.4 35.0 35.0 $95\%$ Cr.I. (38.1, 38.9) (38.1, 38.9)
(38.1, 38.9) (38.1, 38.9) (38.1, 38.9) (35.0, 35.1) (35.0, 35.1) dev† Mean
-1916.3 -1916.3 -1916.3 -1916.3 -1916.3 -1894.0 -1894.0 $95\%$ Cr.I.
(-1920.6,-1908.3) (-1920.6,-1908.3) (-1920.6,-1908.3) (-1920.6,-1908.3)
(-1920.6,-1908.3) (-1899.2,-1885.2) (-1899.2,-1885.2) • ${\dagger}$ deviance
of the model given the data.

The mean and $95\%$ Cr.I. limits of the BMA estimates of parameters of
interest calculated with the use of information criteria score weights
assuming the Gaussian and Inverse Gamma distributions respectively for the
Tetranychus urticae data. Gaussian distribution aic_w dic_w loocv_w waic_w
bic_w elbo_mf elbo_fr $T_{min}$ Mean 10.6 10.6 10.6 10.6 10.6 10.6 10.6 $95\%$
Cr.I. (9.7, 11.4) (9.7, 11.4) (9.7, 11.4) (9.7, 11.4) (9.7, 11.4) (9.7, 11.4)
(9.7, 11.4) $T_{opt}$ Mean 32.6 32.5 32.6 32.6 32.6 32.6 32.6 $95\%$ Cr.I.
(32.0, 33.5) (31.9, 33.4) (32.0, 33.5) (32.0, 33.5) (32.0, 33.5) (32.0, 33.5)
(32.0, 33.5) $T_{max}$ Mean 35.0 35.0 35.0 35.0 35.0 35.0 35.0 $95\%$ Cr.I.
(34.98, 35.0) (34.98, 35.0) (34.98, 35.0) (34.98, 35.0) (34.98, 35.0) (34.98,
35.0) (34.98, 35.0) dev† Mean -818.8 -818.5 -818.7 -818.7 -818.8 -818.8 -818.8
$95\%$ Cr.I. (-833.1, -800.8) (-832.4, -801.1) (-832.9, -800.9) (-832.9,
-800.9) (-833.1, -800.8) (-833.1, -800.8) (-833.1, -800.8) Inverse Gamma
distribution $T_{min}$ Mean 5.0 5.0 5.0 5.0 5.0 11.1 5.0 $95\%$ Cr.I. (4.0,
7.0) (4.0, 7.0) (4.0, 7.0) (4.0, 7.0) (4.0, 7.0) (10.3, 11.9) (4.0, 7.0)
$T_{opt}$ Mean 33.5 33.5 33.5 33.5 33.5 29.3 33.5 $95\%$ Cr.I. (32.3, 34.9)
(32.3, 34.9) (32.3, 34.9) (32.3, 34.9) (32.3, 34.9) (29.2, 29.5) (32.3, 34.9)
$T_{max}$ Mean 33.6 33.6 33.6 33.6 33.6 35.0 33.6 $95\%$ Cr.I. (32.5, 34.9)
(32.5, 34.9) (32.5, 34.9) (32.5, 34.9) (32.5, 34.9) (34.8, 35.0) (32.5, 34.9)
dev† Mean -736.2 -736.2 -736.2 -736.2 -736.2 -563.8 -736.2 $95\%$ Cr.I.
(-742.1, -729.1) (-742.1, -729.1) (-742.0, -729.1) (-742.0, -729.1) (-742.1,
-729.1) (-568.0, -555.9) (-742.1, -729.1) • ${\dagger}$ deviance of the model
given the data.

Posterior summaries for the four models using the Gaussian distribution for
the Propylea Coccinellidae data. In each column we report the HMC, the ADVI-
Mean field and ADVI-Full rank estimates respectively. $T_{min}$ neff‡
$T_{opt}$ neff‡ Mean 95% Cr. I. Bieri 10.6 (9.7, 11.4) 10.7 (10.6, 10.8) 10.7
(10, 11.4) 7906 32.6 (32, 33.5) 32.2 (32, 32.5) 32.2 (31.8, 32.5) 7620 Briere
13.1 (12.4, 13.8) 13.3 (12.9, 13.8) 8.2 (1.8, 15.6) 4 29.7 (29.6, 29.8) 29.7
(29.5, 29.8) 27.3 (11.7, 30.1) 4 Analytis 6.3 (4.1, 10.6) 21.5 (7.7, 55.2) 6.9
(5.4, 9.3) 3742 33.2 (32.1, 34.8) 39.2 (18.1, 89.1) 33.2 (31.7, 37.7) 4435
Lactin -154.4 (-234.7, -119.3) -146.2 (-154.3, -138.1) -159.2 (-243.1, -119.0)
7560 31.2 (30.9, 31.3) 31.1 (31.1, 31.2) 31.1 (30.9, 31.3) 8024 $T_{max}$
neff‡ dev† neff‡ Mean 95% Cr. I. Bieri 35.0 (34.98, 35.02) 33.8 (33.6, 33.9)
33.7 (33.5, 33.9) 7194 -818.7 (-833.1, -800.6) -672.7 (-685.1, -657.4) -672.5
(-686.2, -657.1) 11673 Briere 35.0 (35.0, 35.01) 35.0 (34.8, 35) 32.9 (14.3,
35) 113 -660.5 (-667.0, -649.7) -535.8 (-543.2, -520.5) -429.7 (-535.0,
-402.1) 14478 Analytis 33.6 (32.5, 34.9) 42.5 (32.7, 94.5) 36.1 (4.6, 133.2)
9370 -674.1 (-690.2, -666.3) -160.0 (-261.3, -98.2) -491.9 (-595.9, -293.6)
5164 Lactin 35.0 (35, 35.04) 35.0 (35.0, 35.1) 35.0 (35.0, 35.1) 21622 -809.8
(-823.7, -792.2) -804.7 (-820.5, -784.8) -778.5 (-807.5, -702.6) 12147 •
${\dagger}$ deviance of the model given the data, • ${\ddagger}$ effective
sample size.

Posterior summaries for the four models using the Gaussian distribution for
the Tetranychus urticae data. In each column we report the HMC, the ADVI-Mean
field and ADVI-Full rank estimates respectively. $T_{min}$ neff‡ $T_{opt}$
neff‡ Bieri 10.5 (10.2, 10.8) 9.8 (9.5, 10.2) 9.8 (9.2, 10.3) 15089 33.0
(32.7, 33.4) 158.6 (144.4, 165.4) 183.8 (76.4, 368.5) 8164 Briere 9.3 (8.6,
9.9) 9.3 (9.1, 9.5) 9.3 (8.6, 9.9) 8690 33.1 (32.5, 33.7) 33.0 (32.8, 33.2)
33.1 (32.4, 33.7) 8028 Analytis 4.4 (4.0, 5.3) 4.3 (4.2, 4.4) 5.3 (4.8, 6.1)
8757 32.9 (32.7, 33.3) 33.0 (32.9, 33.1) 33.2 (32.6, 33.7) 41 Lactin 10.4
(10.1, 10.7) -3.9 (-5.8, -1.9) -0.7 (-4.0, 3.0) 11526 32.6 (32.4, 32.8) 32.2
(32.0, 32.3) 32.0 (31.8, 32.2) 9884 $T_{max}$ neff‡ dev† neff‡ Bieri 36.3
(35.9, 36.8) 164.9 (161, 169) 190.6 (80, 376.1) 7561 -1677.1 (-1683.2,
-1666.6) -1461.0 (-1467.7, -1443.6) -1461.7 (-1467.7, -1456.3) 10087 Briere
40.0 (39.3, 40.9) 39.9 (39.6, 40.1) 40.0 (39.1, 40.9) 7562 -1602.2 (-1606.8,
-1593.4) -1593.2 (-1605.7, -1569.1) -1598.5 (-1606.3, -1583.6) 9618 Analytis
35.3 (35.1, 35.5) 35.2 (35.2, 35.3) 35.2 (35.1, 35.4) 8246 -1750.3 (-1758.0,
-1738.3) -1744.7 (-1755.6, -1724) -1709.1 (-1741.9, -1606.9) 10954 Lactin 36.9
(36.6, 37.3) 38.3 (38.1, 38.5) 38.2 (37.9, 38.5) 7797 -1693.9 (-1703, -1679.6)
-1764.0 (-1780.3, -1738.1) -1776.2 (-1784.9, -1762) 8227 • ${\dagger}$
deviance of the model given the data, • ${\ddagger}$ effective sample size.

|  | AIC | DIC | LooCV | WAIC | BIC |  $log\left(P_{y}IS\right)^{\dagger}$ (se) |  $log\left(P_{y}PP\right)^{\ddagger}$ (se) |  $log\left(P_{y}BS\right)$§ (se)
---|---|---|---|---|---|---|---|---|---
Gaussian model | Bieri | -804.7 | -783.6 | -816.1 | -816.1 | -786.1 | 419.9 (19.6) | 433.4 (18.0) | 462.4 (9.2)
Briere | -650.5 | -649.5 | -658.3 | -658.3 | -637.3 | 299.2 (24.6) | 235.4 (6.7) | 231.4 (1.2)
Analytis | -662.1 | -637.1 | -671.4 | -671.4 | -646.1 | 283.7 (24.0) | 312.2 (10.0) | 274.2 (33.7)
Lactin | -789.8 | -776.8 | -806.7 | -806.7 | -763.3 | 390.1 (33.5) | 328.1 (23.0) | 334.7 (1.7)
Inverse Gamma model | Bieri | -716.1 | -719.3 | -719.3 | -719.3 | -702.8 | 334.5 (25.9) | 342.9 (31.5) | 333.2 (5.0)
Briere | -560.8 | -561.4 | -561.4 | -561.4 | -550.2 | 268.2 (23.2) | 265.7 (4.6) | 251.7 (1.9)
Analytis | -734.7 | -732.4 | -732.4 | -732.4 | -718.8 | 340.8 (26.7) | 371.3 (14.5) | 335 (31.2)
Lactin | -698.9 | -699.1 | -699.1 | -699.1 | -685.6 | 329.2 (19.5) | 329.3 (11.8) | 313.7 (18.8)

  * •

${{\dagger}}$ $log\left(P_{y}IS\right)$ denotes the logarithm of estimated
marginal likelihood via Importance sampling,

  * •

${\ddagger}$ $log\left(P_{y}PP\right)$ denotes the logarithm of estimated
marginal likelihood via Power posterior,

  * •

$\S$ $log\left(P_{y}BS\right)$§denotes the logarithm of estimated marginal
likelihood via Bridge sampling.

Table A1: Model selection criteria for the eight models applied to the
Propylea Coccinellidae data

##  Appendix B Posterior prediction graphs for both datasets

Figure 7: Predicted posteriors versus Tetranychus urticae data using Gaussian
distribution Figure 8: Predicted posteriors versus Propylea Coccinellidae data
using Gaussian distribution

Figure 9: Predicted posteriors versus Propylea Coccinellidae data using Zero
Inflated Inverse Gamma distribution

##  Appendix C Bayesian inference

In the current section, we provide details of statistical model specifications
in the Bayesian framework. Specifically, we provide a brief description of the
most important features of Stan’s implementation of HMC and VB so the reader
can get familiar with the tools that Stan is based on. We then provide model
selection and model averaging techniques in order to compare the different
ecological models and to explore and interpret the parameters of interest
combining predictions from all the four of them.

###  C.1 HMC and VB techniques

The HMC method is a Monte Carlo technique that uses Hamiltonian dynamics in
order not only to explore efficiently the target distribution but also to
propose distant samples in the parameter space that do not exclusively depend
on the current state of the Markov chain like considered in previous MCMC
methodology (Neal, 2010). In this way, many performance challenges are tackled
like either the low convergence due to the fact that the parameter space with
high posterior support is not reached or the poor exploration of the target
distribution due to its multi-modality or its shape irregularities. The
existence of Hamiltonian dynamics in the system of the joint density mass
function allows the preservation of volume and hence adequate trajectories can
be used to define complex mappings of the parameter state space without the
need to account for cumbersome Jacobian calculations (Barber et al., 2003).
Thus, by carefully designing automated trajectory realizations in the
Hamiltonian dynamics system, the Stan team managed to create an augmented
software called STAN (Carpenter et al., 2017), which materializes HMC sampling
for the parameters of interest.

Moreover, independently, Automatic Differentiation Variational Inference
(ADVI) technique is referred to the machine learning field (Blei et al.,
2017). The latter is a VB method and posterior target distributions are
approximated by choosing the closest distribution to a parametric family of
tractable distributions like the exponential family via optimization. In order
to achieve this point-wise estimations of the parameters of the family
distribution are estimated so that the Kullback–Leibler ‘KL’ divergence
function is minimized. Specifically, since the KL divergence is intractable
the Evidence Lower Bound is maximized instead (Blei et al., 2017).

###  C.2 Model selection and model averaging

In case there are $m$ models $(\mathcal{M}_{1},\ldots,\mathcal{M}_{m})$ under
consideration, the posterior probability of the suitability of the $i^{th}$
model given the data ${y}$, is given by

$p(\mathcal{M}_{i}|{y})=\frac{p({y}|\mathcal{M}_{i})p(\mathcal{M}_{i})}{\sum_{k=1}^{m}p({y}|\mathcal{M}_{k})p(\mathcal{M}_{k})}$
(10)

where $p(\mathcal{M}_{i})$ expresses the prior belief for the $i^{th}$ model,
while the $p({y}|\mathcal{M}_{i})$ is the model evidence also called ‘marginal
likelihood’ and it can be interpreted as the likelihood over the space of
models, marginalizing out the parameters of the $i^{th}$ model. The ratio of
the marginal likelihoods between two models
$\frac{p({y}|\mathcal{M}_{i})}{p({y}|\mathcal{M}_{j})}$ is called Bayes factor
and is the posterior odds of the null hypothesis that the $i_{th}$ model fits
better the data than the $j_{th}$ model does when the prior probability of the
null is one-half (Kass and Raftery, 1995). The Bayes factor is used in order
to give evidence for the most probable model given the data, when comparing
two alternative models (Kass and Raftery, 1995).

In the case of Bayesian model averaging, the model selection uncertainty is
taking into account in statistical inference. The joint posterior
$p(\mathcal{M}_{i},{\theta}_{i}|{y})$ of the $i^{th}$ model with vector of
parameters $\theta_{i}$, using the Baye’s rule, is proportional to the product
of the likelihood of the $i^{th}$ model times the prior distribution of the
parameters $p({\theta})$ times the prior distribution $p(\mathcal{M}_{i})$
(that expresses our uncertainty of the $i^{th}$ model)

$p(\mathcal{M}_{i},{\theta_{i}}|{y})\propto
p({y}|{\theta}_{i},\mathcal{M}_{i})\cdot p({\theta})\cdot p(\mathcal{M}_{i}).$
(11)

The uncertainty of the $i^{th}$ model given the data can then be re-expressed
via the posterior probability $p(\mathcal{M}_{i}|{y})$ defined by the ratio in
(10), in case of existence of multiple models. The posteriors of the models
can be thought as weights that are critical to the Bayesian model averaging as
they can be used to extract useful weighted statistics from the data
distribution while at the same time taking into account model uncertainties.
Estimation of model parameters and model uncertainties can be achieved either
by directly sampling from the joint posterior (11) or by approximating the
marginal likelihood of each model independently and, accordingly, by
controlling the outcomes with a view to formulating proper weights and
proceeding with the calculation of the averaged statistics. For the former
case, techniques like the reversible jump MCMC (Green, 1995; George and
McCulloch, 1997) and variable selection samples (Carlin and Chib, 1995; Kuo
and Mallick, 1998; Dellaportas et al., 2000) are used. On the other hand, for
the later case, techniques of marginal likelihood approximations via
thermodynamic integration (Friel and Pettitt, 2008), bridge sampling (Meng and
Wong, 1996), importance sampling (Perrakis et al., 2014) or via information-
criteria perspective like in (Kass and Raftery, 1995) are used.

####  C.2.1 Information criteria

The criteria used in the current work are the Akaike information criterion
‘AIC’, the Bayesian information criterion ‘BIC’, the Deviance information
criterion ‘DIC’, the Watanabe–Akaike information criterion ‘WAIC’ and the
Leave-one-out cross-validation criterion ‘Loocv’. Briefly, these criteria
provide an approximation of the expected log predictive density for new-coming
data while correcting bias from data usage. In particular AIC (Akaike, 1974)
is defined as the difference

$AIC(\mathcal{M}_{i})=-2\log(y|\widehat{\theta_{i}})+2k_{i}$

where $\widehat{\theta_{i}}$ is the Maximum Likelihood estimate ‘MLE’ of the
$k_{i}$ parameters of the $i^{th}$ model. Similarly, BIC (Schwarz et al.,
1978) is defined as the difference

$BIC(\mathcal{M}_{i})=-2\log(y|\widehat{\theta_{i}})+k_{i}\cdot\log(n)$

where n is the sample size. In addition, DIC (Spiegelhalter et al., 2002) is
defined as the following difference:

$DIC(\mathcal{M}_{i})=-2\log p(y|\widehat{\theta_{i}})+2p_{DIC}$

where $\widehat{\theta_{i}}$ is the posterior mean of the parameters of the
$i^{th}$ model, whereas $p_{DIC}$ is the effective number of parameters and it
is evaluated following (Spiegelhalter et al., 2002; Gelman et al., 2013) by
either

$p_{DIC_{1}}=E_{\theta|y}\\{-2\log{p}\left(y|\theta\right)\\}+2\log\left\\{{p}\left(y|\widehat{\theta_{i}}\right)\right\\},$

or

$p_{DIC_{2}}=\frac{Var_{\theta|y}\left\\{\log p(y|\theta)\right\\}}{2},$

where $E_{\theta|y}\left\\{\log p(y|\theta)\right\\}$ is an expectation over
the posterior density of $\theta$, whereas $Var_{\theta|y}\left\\{\log
p(y|\widehat{\theta})\right\\}$ is the variance of the log posterior density
of the observed data y, over the posterior density of $\theta$. Furthermore,
WAIC (Watanabe, 2010) is defined as the following difference:

$WAIC(\mathcal{M}_{i})=-2\sum_{j=1}^{N}\log
E_{\theta_{i}}\left\\{p(y_{j}|\theta_{i})\right\\}+2p_{WAIC}$

where $E_{\theta_{i}}\left\\{p(y_{j}|\theta_{i})\right\\}$ is the expectation
of the probability at $y_{j}$ data point over the posterior distribution of
the parameters of the $i^{th}$ model, whereas $p_{WAIC}$ is the effective
number of parameters and it is evaluated following (Gelman et al., 2013) by
either

$p_{{WAIC}_{1}}=2\sum_{j=1}^{N}\left[\log
E_{p(\theta|y)}\left\\{p(y_{j}|\theta)\right\\}-E_{p(\theta|y)}\left\\{\log
p(y_{j}|\theta)\right\\}\right]$

or

$p_{{WAIC}_{2}}=\sum_{j=1}^{N}var_{p(\theta|y)}\left\\{\log
p(y_{j}|\theta)\right\\}$

where $E_{p(\theta|y)}\left\\{\log p(y_{j}|\theta)\right\\}$ is the
expectation over the logarithm of the posterior density of $\theta$ at $y_{j}$
data point, whereas $var_{p(\theta|y)}\left\\{\log p(y_{j}|\theta)\right\\}$
is the variance of the log posterior density of the observed data $y_{j}$,
over the posterior density of $\theta$.

Furthermore, LooCV (Gelman et al., 2013) is defined as the following
difference:

$LooCV(\mathcal{M}_{i})=-2\sum_{j=1}^{N}\log
E_{\theta_{i}^{-j}}\left\\{p(y_{j}|\theta_{i}^{-j})\right\\}-2\beta_{LooCV}$

where $E_{\theta_{i}^{-j}}\left\\{p(y_{j}|\theta_{i}^{-j})\right\\}$ is the
expectation of the probability at $y_{j}$ data point over the posterior
distribution of the parameters of the $i^{th}$ model. The posterior
distribution $p({\theta_{i}^{-j}}|y_{-j})$ is sampled considering a partition
of the data, leaving one data value ($y_{j}$) out of the original sample. The
$\beta_{LooCV}$ is a bias correction of the measure and it is evaluated
following (Gelman et al., 2013) by

$\beta_{Loocv}=\sum_{j=1}^{N}\log
E_{\theta_{i}}\left\\{p(y_{j}|\theta_{i})\right\\}-\frac{1}{N}\sum_{\kappa=1}^{N}\sum_{j=1}^{N}\log
E_{\theta_{i}^{-\kappa}}\left\\{p(y_{j}|\theta_{i}^{-\kappa})\right\\}$

where $E_{\theta_{i}^{-\kappa}}\left\\{p(y_{j}|\theta_{i}^{-\kappa})\right\\}$
is the expectation of the probability at $y_{j}$ data point over the posterior
distribution of the parameters of the $i^{th}$ model leaving out the
$\kappa^{th}$ observation.

####  C.2.2 Marginal likelihood estimation techniques

The marginal likelihood can be viewed as a normalizing constant
$z_{i}=p({y}|\mathcal{M}_{i})$ of the density
$q({\theta_{i}}|{y})=p({y}|{\theta_{i}})\cdot p({\theta_{i}})$ within the
$i^{th}$ ecological model that includes parameters $\theta_{i}$. In the
general scheme of comparing the two densities $q_{0}$ and $q_{1}$ of interest,
as in the case of the Bayes factor of two models or in the case model’s prior
and posterior, a general path from $q_{0}$ to $q_{1}$ can be created according
to (Gelman and Meng, 1998) using a class of densities $p(\theta_{i}|{y},t)$ on
the same space indexed by the continuous auxiliary variable say $t\in[0,1]$. A
key formula that links the corresponding normalizing constant $z(t)$ and the
unormalized density $q(\theta_{i}|{y},t)$ that correspond to the sampling
distribution $p(\theta_{i}|{y},t)$ is given by:

$\frac{d}{dt}\log
z(t)=\int{\frac{1}{z(t)}\frac{d}{dt}q({\theta_{i}}|{y},t)p({\theta_{i}}|t)d{\theta_{i}}}=E_{t}\left\\{\frac{d}{dt}\log{q({\theta_{i}}|{y},t)}\right\\},$
(12)

where the expectation is with respect to the sampling distribution
$p(\theta_{i}|{y},t)$.

In addition, another key formula of estimating a ratio of normalizing
constants has been of great interest such as in computing likelihood ratios in
hypothesis testing or in computational physics in estimating free energy
differences, or in computing the Bayes factor in Bayesian framework (Meng and
Wong, 1996). The general formula is as follows:

$\frac{z_{1}}{z_{0}}=\frac{p({y}|\theta,t=1)}{p({y}|\theta,t=0)}=\frac{E_{0}\left\\{h(\theta)\cdot
q(\theta|{y},t=1)\right\\}}{E_{1}\left\\{h(\theta)\cdot
q(\theta|{y},t=0)\right\\}}$ (13)

where $E_{0}$ and $E_{1}$ expectations are with respect to posterior
distribution densities $p(\theta|{y},t=0)$ and $p(\theta|{y},t=1)$
respectively, whereas the bridge function $h(\theta)$ is defined and
overlapped by the common support of the former densities.

Using general formulas (12) and (13), several marginal probability evaluation
schemes of the $i^{th}$ model are derived (Gelman and Meng, 1998). The power
posterior sampling (Friel and Pettitt, 2008), the importance sampling
(Perrakis et al., 2014) and the bridge sampling (Meng and Wong, 1996;
Overstall and Forster, 2010) techniques are used for the current work.

In the power posterior case, formula (12) is integrated with respect to
variable t and $q(\theta_{i}|{y},t)$ is substituted with density
$p({\theta_{i}}|y)^{t}p(\theta_{i})$. The marginal likelihood
$z_{i}=p({y}|\mathcal{M}_{i})$ is derived from logarithmic scale by the
equation:

$\log\left\\{p({y}|\mathcal{M}_{i})\right\\}=\int_{0}^{1}E_{\theta_{i}|y,t}\left\\{\log{p({y}|\theta_{i})}\right\\}dt$
(14)

where expectation $E_{\theta_{i}|y,t}$ is taken with respect to the density
$p({\theta_{i}}|y)^{t}p(\theta_{i})$ which is defined as the power posterior
at temperature t (Friel and Pettitt, 2008).

Additionally, the standard error ${se}_{i}$ for the ith model estimator (14),
as shown in section E of the appendix is approximated by:

$\hat{{se}_{i}}=\sqrt{\frac{\left(t_{2}-t_{1}\right)^{2}}{2}s_{1}^{2}+\sum_{k=2}^{n-1}\frac{\left(t_{k}-t_{k-1}\right)^{2}}{2}s_{k}^{2}+\frac{\left(t_{n}-t_{n-1}\right)^{2}}{2}s_{n}^{2}},$

where $t_{k}$ is the time after discretization $0=t_{0}<t_{1}<t_{k}<t_{n}=1$
and $s_{k}$ is the standard error of the corresponding estimation
$\log\left\\{p({y}|\mathcal{M}_{i})\right\\}$ given in (14).

In the case of importance sampling, the marginal likelihood is assessed by
introducing the proper density function g. After sampling from the proposed
density function g, the marginal likelihood is calculated as with respect to g
as:

$p({y}|\mathcal{M}_{i})=E_{g}\left\\{\frac{q(\theta_{i}|{y})}{g(\theta_{i})}\right\\}$

Following (Perrakis et al., 2014), we use the density $q(\theta_{i}|{y})$
equal to $p({y}|\theta_{i},\phi_{i})\cdot p(\theta_{i},\phi_{i})$ and the
auxiliary importance function g used is as follows:

$g(\theta_{i})=g(\theta_{i},\phi_{i})=p(\theta_{i}|{y})p(\phi_{i}|{y}),$ (15)

where $(\theta_{i},\phi_{i})$ are the parameters of the $i^{th}$ model divided
into two blocks $\theta_{i}$ and $\phi_{i}$ which may or may not be
independent. The right hand side of (15) is the product of the marginal
posterior distributions of the block. Thus, the marginal probability which
gives the target value is given as follows:

$p({y}|\mathcal{M}_{i})=\iint\frac{p({y},\theta_{i},\phi_{i})}{g(\theta_{i},\phi_{i})}g(\theta_{i},\phi_{i})d(\theta_{i},\phi_{i})=E_{g}\left\\{\frac{p({y},\theta_{i},\phi_{i})}{g(\theta_{i},\phi_{i})}\right\\}$
(16)

The standard error ${se}_{i}$ of (16) as shown in section F of the appendix
is:

${\hat{{se}_{i}}}=\sqrt{\frac{1}{K}\sum_{j=1}^{K}\left\\{\frac{p({y}|\theta^{j},\phi^{j})\cdot
p(\theta^{j},\phi^{j})}{g\left(\theta^{j}\right)}-{\hat{z}}_{i}\right\\}^{2}},$

where $\widehat{z}_{i}$ is the estimation of the corresponding marginal
probability (in the same form of (16)), while $(\theta^{j},\phi^{j})$ are
draws $j=1,2,\dots,K$ from the importance function in (15).

Additionally, using an alternative version of (13) in (Meng and Wong, 1996;
Frühwirth-Schnatter, 2004; Overstall and Forster, 2010) the marginal
likelihood of a single model is evaluated using bridge sampling by the
formula:

$z_{i}=p({y}|\mathcal{M}_{i})=\frac{E_{g}\left\\{h(\theta_{i})\cdot
q(\theta_{i}|{y})\right\\}}{E_{p}\left\\{h(\theta_{i})\cdot
g(\theta_{i})\right\\}},$ (17)

where $E_{g}$ and $E_{p}$ are the expectations with respect to $g(\theta_{i})$
a so-called proposal distribution and to $p(\theta_{i}|{y})$ the $i^{th}$
model posterior distribution respectively.

The bridge function $h(\theta_{i})$ is selected to minimize the relative mean-
squared error of (13). Following (Meng and Wong, 1996) the bridge function is
specified by:

$h(\theta_{i})=\mathcal{C}\cdot\frac{1}{s_{1}\cdot
q(\theta_{i}|{y})+s_{2}\cdot p(y)\cdot g(\theta_{i})},$ (18)

where $s_{1}=\frac{N_{1}}{N_{1}+N_{2}}$, $s_{2}=\frac{N_{2}}{N_{1}+N_{2}}$ and
$\mathcal{C}$ is a constant. $N_{1}$ is the sample size from the posterior and
$N_{2}$ is the sample size from $g(\theta_{i})$.

The optimal bridge function in (18) includes the marginal likelihood under-
assessment so that it cannot be evaluated directly. For this purpose the
iterative method suggested by (Meng and Wong, 1996) and applied in (Gronau et
al., 2020) in R software (R Core Team, 2020) is used. The alternatives used in
place of distribution $g$ is either a multivariate normal distribution with
mean vector and covariance matrix that match the respective posterior samples
quantities or a standard multivariate normal distribution in combination with
a warped posterior distribution of which the first three moments correspond to
(Gronau et al., 2020).

Moreover, following (Frühwirth-Schnatter, 2004) the relative mean square error
${RE}^{2}_{i}=\frac{E\left\\{\widehat{z_{i}}-z_{i}\right\\}^{2}}{z^{2}_{i}}$
of (17) is evaluated by the formula:

$\widehat{RE}^{2}_{i}={\frac{1}{N_{2}}\frac{V_{g}\left\\{f_{1}(\theta_{i})\right\\}}{E_{g}^{2}\left\\{f_{1}(\theta_{i})\right\\}}+\frac{\rho_{f_{2}}(0)}{N_{1}}\frac{V_{p}\left\\{f_{2}(\theta_{i})\right\\}}{E_{p}^{2}\left\\{f_{2}(\theta_{i})\right\\}}},$
(19)

where $f_{1}(\theta_{i})=\frac{q(\theta_{i}|{y})}{s_{1}\cdot
q(\theta_{i}|{y})+s_{2}\cdot g(\theta_{i})}$,
$f_{2}(\theta_{i})=\frac{g(\theta_{i})}{s_{1}\cdot
q(\theta_{i}|{y})+s_{2}\cdot g(\theta_{i})}$,
$V_{g}(f_{1}(\theta_{i}))=\int\left\\{f_{1}(\theta_{i})-E(f_{1}(\theta_{i}))\right\\}^{2}g(\theta_{i})d\theta$
is the variance of $f_{1}(\theta_{i})$ with respect to $g(\theta_{i})$. The
term $\rho_{f_{2}}(0)$ in (19) corresponds to the normalized spectral density
of the auto-correlated process $f_{2}(\theta_{i})$ at the frequency 0.

Following (Gronau et al., 2020) the square root of ${RE}^{2}$ can be
interpreted as coefficient of variation provided that the bridge sampling
estimator $\widehat{z}_{i}$ is unbiased. Then the stantard error ${se}_{i}$ of
the bridge estimator is evaluated by the product
$\widehat{se}_{i}=\widehat{RE}\cdot E(\widehat{z}_{i})$

####  C.2.3 BMA weights

We can derive a weighted prediction $\tilde{y}$ over the m different models
$\mathcal{M}_{1},\mathcal{M}_{2},\dots\mathcal{M}_{m}$ predictions
$\hat{y}_{1},\hat{y}_{2},\dots\hat{y}_{m}$ by imposing appropriate weights
$w_{1},w_{2},\dots w_{m}$.

$\tilde{y}=\sum_{i=1}^{m}{\hat{y}_{i}\cdot
w_{i}}\quad\textrm{and}\quad\sum_{i=1}^{m}{w_{i}}=1.$

In the Bayesian framework, model weights definition is straightforward. The
model weights used are the posterior model weights
$w_{i}=p(\mathcal{M}_{i}|{y})$ given in (10) that represent the relative
probability of each model given the data. So a major challenge is to estimate
these Bayesian weights. Except for using the marginal likelihood estimations
mentioned in previous section, we also investigate approximations of the
weights by using the BIC for each model. In particular, model weights can be
estimated through the following equations (Kass and Raftery, 1995; Buckland et
al., 1997):

(20)
$\displaystyle{w}_{i}=\frac{e^{-0.5\cdot\left(BIC(\mathcal{M}_{i})\right)}}{\sum_{j=1}^{m}{e^{-0.5\cdot\left(BIC(\mathcal{M}_{j})\right)}}}.$

Instead of BIC, the AIC, DIC, WAIC and LooCV are also used in (20). We
investigate both approaches in insect observed rates and compare the results
taking into account model complexity, data scarcity and Biological
interpretation.

##  Appendix D Power posterior for Gaussian and Inverse Gamma distribution

Following (Friel and Pettitt, 2008) the power posterior in (14) includes the
likelihood raised to the power of t $p({\theta_{i}}|y)^{t}$.

In the the Gaussian case, the likelihood involved becomes:

$\displaystyle
p(y|{\theta})^{t}=p^{t}\left(y|\mu,\sigma^{2}\right)=\frac{e^{-\frac{1}{2}\left(\sqrt{t}\cdot\frac{y-\mu}{\sigma}\right)^{2}}}{\sqrt{2\pi\sigma^{2}}}=p\left(y|\mu,\frac{\sigma^{2}}{t}\right)\cdot{\frac{\sqrt{t}}{t}}.$

Inserting the current ecological model $r\left(T;\theta\right)$, the log of
the power posterior is given by:

$\displaystyle
log\left(p(y|{\theta})^{t}\right)=-\frac{1}{2}\log{\left(2\pi\sigma^{2}\right)}-\frac{1}{2}\left(\sqrt{t}\cdot\frac{y-r\left(T;\theta\right)}{\sigma}\right)^{2}.$

In the Inverse Gamma case, the likelihood involved becomes:

$\displaystyle
p(y|{\theta})^{t}=p(y|{\alpha,\beta})^{t}=\frac{\beta^{\alpha\cdot
t}}{\Gamma^{t}\left(\alpha\right)}y^{-\alpha\cdot
t-t}exp\left(-\frac{\beta\cdot t}{y}\right)=$
$\displaystyle\frac{\left(\beta\cdot t\right)^{\alpha\cdot
t}}{\Gamma\left(\alpha\cdot t\right)}\frac{\Gamma\left(\alpha\cdot
t\right)}{\Gamma^{t}\left(\alpha\right)\cdot(t)^{\alpha\cdot
t}}y^{-\alpha\cdot t-1}y^{1-t}exp\left(-\frac{\beta\cdot t}{y}\right)=$
$\displaystyle InG\left(\alpha\cdot t,\beta\cdot
t\right)\cdot\frac{\Gamma\left(\alpha\cdot
t\right)}{\Gamma^{t}\left(\alpha\right)\cdot(t)^{\alpha\cdot t}}\cdot
y^{1-t}.$

Inserting the current ecological model $r\left(T;\theta\right)$, the log of
the power posterior is given by:

$\displaystyle log\left(p(y|{\theta})^{t}\right)={\alpha\cdot
t}\cdot\left(log\left(\alpha-1\right)+log\left(r\left(T;\theta\right)\right)\right)-$
$\displaystyle t\cdot
log\left(\Gamma\left(\alpha\right)\right)-t\cdot\left(\alpha+1\right)\cdot
log\left(y\right)-\frac{\left(\alpha-1\right)\cdot\left(r\left(T;\theta\right)\right)\cdot
t}{y}.$

##  Appendix E Estimation of Variance of Power posterior method

$\displaystyle\hat{\sigma_{y}}^{2}=Var\left\\{logP\left(y|\theta\right)\right\\}=Var\left\\{\int_{0}^{1}{E\left\\{{log}P{\left(y|\theta\right)}\right\\}dt}\right\\}\simeq$
$\displaystyle\sum_{i=1}^{n-1}{\frac{1}{4}Var\left\\{E_{\theta|y,t_{i}}\left\\{{log}P{\left(y|\theta\right)}\right\\}+E_{\theta|y,t_{i}+1}\left\\{{log}P{\left(y|\theta\right)}\right\\}\right\\}\left(t_{i+1}-t_{i}\right)^{2}}\simeq$
$\displaystyle\sum_{i=1}^{n-1}{\left[\frac{{{sd}^{2}}_{\theta|y,t_{i}}}{2}\right]\cdot\left(t_{i+1}-t_{i}\right)^{2}}$

where ${{sd}}_{\theta|y,t_{i}}$ is the std error estimated at the $t_{i}$
temperature

##  Appendix F Estimation of Variance of Importance sampling method

The marginal likelihood estimate for the ith model is given by
$m=p\left(y|M_{i}\right)=E_{g}\left\\{\frac{q\left(\theta_{i}|y\right)}{g\left(\theta_{i}\right)}\right\\}=\int\frac{q\left(\theta_{i}|y\right)}{g\left(\theta_{i}\right)}g\left(\theta_{i}\right)d\theta_{i}$
provided $g\left(\theta-i\right)>0$ whenever $q\left(\theta_{i}|y\right)\neq
0$ where the density $q\left(\theta_{i}\middle|\ y\right)$ is equal to
$p\left(y\middle|\theta_{i},\phi_{i}\right)\cdot
p\left(\theta_{i},\phi_{i}\right)$ and the auxiliary importance function $g$
used is the following:
$g\left(\theta_{i}\right)=g\left(\theta_{i},\phi_{i}\right)=p\left(\theta_{i}\middle|\
y\right)p\left(\phi_{i}\middle|\ y\right)$
Removing the ith model index we can evaluate the marginal likelihood via MC
integration which gives the formula below:

${\hat{z}}_{y}=\hat{p}\left(y\right)=\frac{1}{M}\sum_{J=1}^{M}\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}$

where $\theta^{j},j=1,2,3,\ldots,M$ are obtained from density:
$g\left(\theta_{i}\right)$.

$E_{g}\left({\hat{z}}_{y}\right)=\frac{1}{M}\sum_{j=1}^{M}{E_{g}\left\\{\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}\right\\}}=p\left(y\right)=m.$

Assuming that
$Cov\left\\{\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)},\frac{q\left(\theta^{i}|y\right)}{g\left(\theta^{i}\right)}\right\\}=0,\
for\ i\neq\ j$

$V_{g}\left({\hat{z}}_{y}\right)=V_{g}\left(\frac{1}{M}\sum_{j=1}^{M}\left\\{\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}\right\\}\right)=\frac{1}{M^{2}}\sum_{j=1}^{M}{V_{g}\left\\{\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}\right\\}+}\frac{1}{M^{2}}\sum_{k,j=1}^{M}{{Cov}_{g}\left\\{\frac{q\left(\theta^{k}|y\right)}{g\left(\theta^{k}\right)},\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}\right\\}=}$
$\frac{1}{M^{2}}\sum_{j=1}^{M}{V_{g}\left\\{\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}\right\\}}=\frac{1}{M}V_{g}\left\\{\frac{q\left(\theta|y\right)}{g\left(\theta\right)}\right\\}=\frac{1}{M}\left\\{{E_{g}\left\\{\frac{q\left(\theta|y\right)}{g\left(\theta\right)}\right\\}}^{2}-{E_{g}}^{2}\left\\{\frac{q\left(\theta|y\right)}{g\left(\theta\right)}\right\\}\right\\}=$
$\frac{1}{M}\int_{G}{\left\\{\frac{q\left(\theta|y\right)}{g\left(\theta\right)}\right\\}^{2}g\left(\theta\right)d\theta}-\frac{1}{M}m^{2}\Rightarrow$
Also, using properties of expectation and the fact that $g\left(\theta\right)$
is proper
$V_{g}\left({\hat{z}}_{y}\right)=\frac{1}{M}\int_{G}{\left\\{\frac{q\left(\theta|y\right)-m\cdot
g\left(\theta\right)}{g\left(\theta\right)}\right\\}^{2}g\left(\theta\right)d\theta}=\frac{\sigma_{y}^{2}}{M}$
The latest can be used to estimate the std error of the estimator
${\hat{z}}_{y}$ via the formula:

${\hat{\sigma}}_{y}=\sqrt{\frac{1}{K}\sum_{j=1}^{K}\left\\{\frac{q\left(\theta^{j}|y\right)}{g\left(\theta^{j}\right)}-{\hat{z}}_{y}\right\\}^{2}}$

where $\theta^{j},j=1,2,3,\ldots,K$ are obtained from density
$g\left(\theta_{i}\right)$.