# Bayesian inference with tmbstan for a state-space model with VAR(1) state
equation

Yihan Cao Centre for Biodiversity Dynamics, Department of Mathematical
Sciences, Norwegian University of Science and Technology, 7491 Trondheim,
Norway Jarle Tufto11footnotemark: 1

##  1 Introduction

Both frequentist and Bayesian statistical inference have been used for
investigating ecological processes. In the frequentist framework, Template
model builder <TMB,>kristensen2016tmb, an R package developed for fast fitting
complex linear or nonlinear mixed models, has gained the popularity recently,
especially in the field of ecology which usually involves in modeling
complicated ecological processes <for
example>cadigan2015state,albertsen2016choosing,auger2017spatiotemporal. The
combination of reverse-mode automatic differentiation and Laplace
approximation for high-dimension integrals makes parameter estimation with TMB
very efficient even for non-Gaussian and complex hierarchical models. TMB
provides a flexible framework in model formulation and can be implemented even
for statistical models where the predictor is nonlinear in parameters and
random effect. However, the lack of capability of working in the Bayesian
framework has hindered the adoption of it for Bayesians.

Within the Bayesian framework, the software package Stan Gelman . (2015), a
probabilistic programming language for statistical inference written in C++
attracts peoples attention. It uses the No-U-Turn Sampler (NUTS) Hoffman
Gelman (2014), an adaptive extension to Hamiltonian Monte Carlo Neal . (2011),
which itself is a generalization of the familiar Metropolis algorithm, to
conduct sampling more efficiently through the posterior distribution by
performing multiple steps per iteration. Stan is a valuable tool for many
ecologists utilizing Bayesian inference, particularly for problems where BUGS
Lunn . (2000) is prohibitively slow Monnahan . (2017). As such, it can extend
the boundaries of feasible models for applied problems, leading to a better
understanding of ecological processes. Fields that would likely benefit
include estimation of individual and population growth rates, meta-analyses
and cross-system comparisons, among many others.

Combining the merits of TMB and Stan, the new software package tmbstan
Monnahan Kristensen (2018) which provides MCMC sampling for TMB models was
developed. This package provides ADMB and TMB users a possibility for making
Bayesian statistical analysis when prior information on the unknown parameters
is available. From the user’s perspective, it implements NUTS sampling from a
target density proportional to the product of marginal likelihood (computed by
TMB or Stan) and the prior density specified by the user. The user has the
flexibility to decide which random effects are integrated out via the Laplace
approximation in TMB and then the TMB model is passed to function Stan in the
RStan package so that the rest of the parameters are integrated by Stan. This
methodology might therefore potentially be more computationally efficient than
using MCMC alone to integrate out all parameters. monnahan2018no introduced
the tmbstan package, applied it to simulation studies and compared its
capabilities (computational efficiency and the accuracy of Laplace
approximation) with other platforms such as ADMB and TMB.

However, it is unclear that if Bayesian inference with arbitrary prior
distribution implemented with Stan would perform comparatively with
frequentist inference when modeling complex ecological processes. It is also
unclear that when using tmbstan, if using the Laplace approximation to
integrate latent variables is more computationally efficient than handling all
latent variables via MCMC. In the case studies in monnahan2018no, Laplace
approximation turned out to reduce the computational efficiency of MCMC.
Another issue arose in the case studies is that the Laplace approximation to
the integration of random effects is not accurate to a degree and this could
lead to biased parameter estimates or uncertainties in parameter estimation.
To gain more insights on these issues, in this paper we conduct simulation
studies and a case study in the context of modeling fluctuating and auto-
correlated selection with state-space models (SSM). These forms of models are
more generally increasingly used in ecology to model time-series such as
animal movement paths and population dynamics <for
example>cadigan2015state,albertsen2016choosing,auger2017spatiotemporal.
Furthermore, following cao2019time, we also use order-1 vector autoregressive
model (VAR(1)) to model the unobserved states, which in our study are
temporally fluctuating and potentially auto-correlated height, width and
location of a Gaussian fitness function. This also allows us to make a further
investigation into the issue of underestimation of the auto-correlation
parameter in auto-regressive models shown in Chevin2015 and cao2019time.

Through the simulation and empirical studies, our paper aims to (1) compare
estimates between frequentist inference and Bayesian inference under different
simulation schemes; (2) investigate how the choice of prior influence Bayesian
inference; (3) compare the computational efficiency of MCMC with and without
integrating out some of the random effects via Laplace approximation.

##  2 Methodology

###  2.1 Model formulation

We consider a typical ecological process, the fluctuating selection in a bird
species, the great tit (Parus major). We conduct the study in the context of
temporally changing selection on the laying date with the number of fledglings
as the fitness component, but it can be generalized to any episode of
viability or fertility selection, or to overall selection through lifetime
fitness. The discrete nonnegative variable, number of fledglings, is best
modelled by distributions such as Poisson, or zero-inflated Poisson <for
example>Chevin2015,cao2019time. Within the framework of generalized linear
models, the expected value of response variable is commonly linked to the
linear predictors of biologically interest by logarithm. When both linear and
quadratic effects of the traits are included, this leads to a Gaussian model
of stabilizing selection. In this study, the number of fledglings in a
specific brood is assumed to be Poisson distributed,
$X_{i}|w_{i}\sim\operatorname{Poisson}(w_{i})$, where $i$ indicates the
breeding event. The fitness (the expected number of fledglings $w_{i}$ ) of
individuals with phenotype $z_{i}$ is then given by

$\displaystyle\ln w_{i}$
$\displaystyle=\eta_{t}^{(\alpha)}-\frac{(z_{i}-\eta_{t}^{(\theta)})^{2}}{2(e^{\eta_{t}^{(\omega)}})^{2}},$
(1)

where $\eta_{t}^{(\alpha)}$, $\eta_{t}^{(\theta)}$ and
$e^{\eta_{t}^{(\omega)}}$ ($e$ based to guarantee positive) are parameters
determining the logrithm of maximum fitness, optimum laying date and width of
the fitness function in year $t$ respectively. We further model
$\eta_{t}^{(\alpha)}$, $\eta_{t}^{(\theta)}$ and $\eta_{t}^{(\omega)}$, the
three stochastic processes as following:

$\displaystyle\eta_{t}^{(\alpha)}$
$\displaystyle=\mu_{\alpha}+\sigma_{\alpha}\alpha_{t},$ (2)
$\displaystyle\eta_{t}^{(\theta)}$
$\displaystyle=\mu_{\theta}+\sigma_{\theta}\theta_{t},$
$\displaystyle\eta_{t}^{(\omega)}$
$\displaystyle=\mu_{\omega}+\sigma_{\omega}\omega_{t}.$

The elements of vector
$\mathbf{\mu}=(\mu_{\alpha},\mu_{\theta},\mu_{\omega})^{T}$ are the means of
the three processes. The stochastic processes $\alpha_{t}$, $\theta_{t}$,
$\omega_{t}$ are assumed to be multivariate normal distributed
$(\alpha_{t},\theta_{t},\omega_{t})^{T}\sim\mathbf{N}_{3}(\mathbf{0},\mathbf{\Gamma_{0}})$
with
$\mathbf{\Gamma_{0}}=\begin{bmatrix}1&\rho_{\alpha,\theta}&\rho_{\alpha,\omega}\\\
\rho_{\alpha,\theta}&1&\rho_{\theta,\omega}\\\
\rho_{\alpha,\omega}&\rho_{\theta,\omega}&1\end{bmatrix}$, where
$\rho_{\alpha,\theta}$, $\rho_{\alpha,\omega}$ and $\rho_{\theta,\omega}$
indicate the correlations and are assumed to be mutually independent.
$(\alpha_{t},\theta_{t},\omega_{t})^{T}$ are further assumed to follow a
first-order vector autoregressive (VAR(1)) process as below:

$\left[\begin{matrix}\alpha_{t}\\\ \theta_{t}\\\
\omega_{t}\end{matrix}\right]=\mathbf{\Phi}\left[\begin{matrix}\alpha_{t-1}\\\
\theta_{t-1}\\\ \omega_{t-1}\end{matrix}\right]+\mathbf{w}_{t},$ (3)

where $\mathbf{\Phi}$ is $3\times 3$ transition matrix and $\mathbf{w}_{t}$ is
a 3-dimentional vector of white noise. The covariance matrix of
$\mathbf{w}_{t}$ is calculated as
$\mathbf{\Gamma_{0}}-\mathbf{\Phi}\mathbf{\Gamma_{0}}\mathbf{\Phi}$.
Correlations between the elements of $\mathbf{w}_{t}$ are determined by both
$\mathbf{\rho}=(\rho_{\alpha,\theta},\rho_{\alpha,\omega},\rho_{\theta,\omega})$
and $\mathbf{\Phi}$. If $\mathbf{\rho}$ is $\mathbf{0}$ vector and
$\mathbf{\Phi}$ is diagonal, then $\mathbf{w}_{t}$ reduces to be three
independent and identically distributed white noise processes. In this case,
$\alpha_{t}$, $\theta_{t}$ and $\omega_{t}$ simplify to three independent
first-order autoregressive (AR(1)) processes. If $\rho$ is $\mathbf{0}$ and
all entries of $\mathbf{\Phi}$ are zero, both
$(\alpha_{t},\theta_{t},\omega_{t})^{T}$ and $\mathbf{w}_{t}$ reduce to three
independent and identically distributed white noise processes. In any case,
our non-centered parameterization implies that the standard deviation of
$\alpha_{t}$, $\theta_{t}$ and $\omega_{t}$ is only determined by
$\sigma_{\alpha}$, $\sigma_{\theta}$ and $\sigma_{\omega}$ respectively. We
expect the non-centered parameterization yields simpler posterior geometries
Betancourt Girolami (2015) and will be much more efficient in terms of
effective sample size when there is not much data (Stan Development Team,
20182, chapter 20).

It is worth mentioning that one objective of this study is to provide another
case study beyond the ones in monnahan2018no. Therefore, even though
$\alpha_{t}$, $\theta_{t}$ and $\omega_{t}$ are assumed to be VAR(1) in the
model, in the simulation study we consider only AR(1) $\theta_{t}$ and white
noise of $\alpha_{t}$ and $\omega_{t}$. The alternative simulation studies in
which $\alpha_{t}$, $\theta_{t}$ and $\omega_{t}$ are formulated as other
possible stochastic processes can be conducted similarly and exhaustively, but
that is an enormous amount of work in one single study. When $\alpha_{t}$,
$\theta_{t}$ and $\omega_{t}$ are assumed to be VAR(1), one caution to be
taken is that all the eigenvalues of $\mathbf{\Phi}$ must lie in the unit
circle to guarantee the VAR (1) process to be stationary Wei (2006). At last,
in the simulation study, we assume that the model structure is known, which
means that we already know $\theta_{t}$ is AR(1) process since the aim of the
study is not to explore the structure of the true model.

###  2.2 Prior distribution

The priors are assumed to be independent to each other
$\pi(\mathbf{\mu,\Phi,\Sigma})=\pi(\mathbf{\mu})\pi(\mathbf{\Phi})\pi(\mathbf{\Sigma})$.
We take a normal $N(\mathbf{m},q\mathbf{I}_{3})$ prior distribution for the
process mean vector $\mathbf{\mu}=(\mu_{\alpha},\mu_{\theta},\mu_{\omega})$
and input weak prior information on the process mean by taking
$\mathbf{m}=\mathbf{0}$ and $q=100$. Since in this study we assume constant
$\eta_{t}^{(\alpha)}$ and $\eta_{t}^{(\omega)}$, $\phi_{\theta,\theta}$ is the
only non-zero entry in $\mathbf{\Phi}$. We used truncated normal prior on
$\phi_{\theta,\theta}$ since it outperforms Jeffreys’ prior Jeffreys Jeffreys
(1961), g prior Zellner (1986) and natural conjugate prior Schlaifer Raiffa
(1961) in terms of posterior sensitivity using Highest Posterior Density
Region (HPDR) criterion concluded from the simulation study in
karakani2016bayesian. lei2011bayesian also uses truncated normal distribution
as subjective prior for the auto-regressive parameter in its AR (1) model. The
mean and standard deviation of the truncated normal distribution are
arbitrarily set to be 0 and 0.5 respectively.

For the variance of the error term $\sigma_{\theta}^{2}$
($\sigma_{\alpha}^{2}$ and $\sigma_{\omega}^{2}$ are assumed to be zero), two
priors are used:
(1) half-Cauchy (0, 10) prior on $\sigma_{\theta}$ (Prior1);
(2) lognormal (1, 0.5) prior on $\sigma_{\theta}$ (Prior2).
These two priors are referred to Prior1 and Prior2 respectively in the rest of
this paper. It is worth mentioning that we also tested uniform prior on
$log(\sigma_{\theta})$ (non-informative improper prior which equals to
$1/\sigma$ prior on $\sigma$ Gelman . (2006)) and inverse-gamma (1, 1) prior
on $\sigma_{\theta}^{2}$ (non-informative proper prior, also illustrated in
Gelman . (2006)), but both of them render an issue that the sampler traps in a
subspace of the whole parameter space of $log(\sigma_{\theta})$ and results in
numerous divergent transitions. It was potentially caused by the posterior
becoming improper and consisting of a mode and an infinite low-posterior-
density ridge extending to infinity as illustrated in tufto2012estimating. We
thus in this study only consider the two proper informative priors (Prior1 and
Prior2), while more information on the MCMC with inverse-gamma (1, 1) prior on
$\sigma_{\theta}^{2}$ is given in Supporting Information.

Note also that the scale parameters $\log(\sigma_{\theta})$ is declared in the
TMB template in the logarithmic format, but the half-Cauchy prior and
lognormal prior contributed to the total likelihood with the log density in
terms of $\sigma_{\theta}$ and for inverse-gamma prior, it is in terms of
$\sigma_{\theta}^{2}$, where $\sigma_{\theta}$ is a positive transform
$\sigma=e^{log\sigma}$. Therefore, Jacobian adjustment (see chapter 20.3 in
stan2018stan for Jacobian adjustment) was conducted by adding
$log\sigma_{\theta}$ to the total likelihood when half-Cauchy prior and
lognormal prior are used. When testing inverse-gamma prior, it was $\log
2+2\log\sigma_{\theta}$ added to the total likelihood.

###  2.3 Software implementation

The model is formulated with C++ and passed to TMB for frequentist inference.
The model objective (fn) and gradient (gr) functions are fed to optimization
function nlminb with default setting to optimize the objective function.

For Bayesian inference, the TMB model objective and gradient functions are
passed to tmbstan which uses the stan function and executes the No-U-Turn
sampler (NUTS) algorithm by default to sample. Currently the other options are
"HMC" (Hamiltonian Monte Carlo), and "Fixed param". We ran the simulation
study on a multicore computing server with enough RAM to avoid swapping to
disk. The number of warmup iterations to be excluded when computing the
summaries is set to 1000 and for total sample length, it is 3000. We thin each
chain to every second sample and set the value adapt delta to 0.95, which is
the average proposal acceptance probability Stan aims for during the adaption
(warmup) period. We set a seed for each simulation including data set and
tmbstan to make sure all the simulation results are reproducible.

Divergent transitions during sampling may occur due to a large step size in
the sampler or a poorly parameterized model, meaning that the iteration of the
MCMC sampler runs into numerical instabilities Carpenter . (2017) and thus
inferences will be biased. RStan team suggested that the problem may be
alleviated by increasing the adapt delta parameter (gives a smaller step
size), especially when the number of divergent transitions is small Stan
Development Team (20181). In our simulation studies, we find it difficult to
completely avoid divergent transitions across all data sets even though adapt
delta is increased to 0.95. Similar to fuglstad2019intuitive, we thus removed
simulations where 0.1% or more divergent transitions in the iterations after
warmup occur during the inference to avoid reporting biased results.

It is worth mentioning that the execution of Markov chains can be done in
parallel. While the default of RStan is to use 1 core, the RStan team
recommended to set it to as many processors as the hardware and RAM allow and
at most one core per chain Stan Development Team (20181). The simulations we
run are done with a server that has 28 available cores. We thus set the number
of cores to be 4 for the 4 Markov chains. However, since R is single threaded
and for frequentist inference, optimization algorithm used in R function
"nlminb" only uses one core of CPU, we thus only compare the computational
efficiency between tmbstan with and without Laplace approximation and ignore
the computational efficiency with "nlminb" to ensure fair comparisons.

##  3 Simulation scheme and results

###  3.1 Simulation scheme

All the data simulated are in natural units and considered to be biologically
realistic according to the empirical studies of natural birds populations
<e.g.>grant2002unpredictable, vedder2013quantitative. Samples were modeled
from a population undergoing stabilizing selection with AR(1) $\theta_{t}$,
fixed $\eta_{t}^{(\alpha)}$ and $\eta_{t}^{(\omega)}$. Vector
$\mathbf{\mu}=(\mu_{\alpha},\mu_{\theta},\mu_{\omega})^{T}$ is set to
$(2,20,3.5)$. The autocorrelation $\phi_{\theta,\theta}$ is set to 0.1, 0.4
and 0.7 (only positive values considered since the estimate of auto-
correlation in temporal optimal laying date is positive, for example 0.3029 in
Chevin2015 and 0.524 in cao2019time), the variance of fluctuating optimal
laying date $\sigma_{\theta}$ is set to 20.

For each value of $\phi_{\theta,\theta}$, $tmax=25$ or 50 time points were
simulated and for each time point the sample size was drawn from a Poisson
distribution with mean $n=25$, 50 or 100 individuals. We considered four
combinations of $tmax$ and $n$, which are $(tmax=25,n=50)$, $(tmax=25,n=100)$,
$(tmax=50,n=25)$ and $(tmax=50,n=100)$. These four combinations are refered as
simulation setting 1, 2, 3, 4 respectively in the following sections. Similar
to cao2019time, we neglected response to selection and used the same normal
distribution for simulating individual phenotype each year. The phenotypic
standard deviation before selection $\sigma_{z}$ was set to 20, such that the
strength of stabilizing selection
$S=\sigma^{2}_{z}/(e^{\eta_{t}^{(\omega)}})^{2}+\sigma^{2}_{z}$
<e.g.>Chevin2015 was 0.267. For each individual, its fitness was computed from
its phenotype using equation (1), and its actual number of offspring was then
drawn from a Poisson distribution with mean $w_{t}(z)$.

###  3.2 Frequentist vs. Bayesian estimates

The results of one single simulation obtained from maximum likelihood in the
frequentist framework are compared with those from tmbstan. The summaries of
the estimates with tmbstan are computed after dropping the warmup iterations
and merging the draws from all the four chains. The frequentist and Bayesian
estimates with different sample sizes and $\phi_{\theta,\theta}=0.4$ are shown
in Table 1, the estimates with other values of auto-correlation in
$\theta_{t}$ ($\phi_{\theta,\theta}$ =0.1 and 0.7) can be found in Supporting
Information.

Table 1: Frequentsit and Bayesian estimates (standard errors) from the model
with AR(1) $\theta_{t}$, autocorrelation in $\theta_{t}$
$\phi_{\theta,\theta}=0.4$, and different sample sizes ($(tmax=25,n=50)$,
$(tmax=25,n=100)$, $(tmax=50,n=25)$ and $(tmax=50,n=100)$) from one
realization of the simulation. For each sample size setting, the number of
divergent transitions in the MCMC is also reported and is used as a measure of
stability of the inference scheme. MLE stands for maximum likelihood estimate,
Prior1 and Prior2 represent half-Cauchy (0, 10) and lognormal (1, 0.5) prior
respectively.  $\phi_{\theta,\theta}=0.4$, $tmax=25$, $n=50$
---
Parameters | True value | MLE | Prior1 | Prior2
no. divergent transitions | NA | NA | 1 | 0
$\mu_{\alpha}$ | 2 | 2.017(0.015) | 2.017(0.015) | 2.016(0.015)
$\mu_{\theta}$ | 20 | 18.5(3.7) | 18.3(5.1) | 18.5(3.7)
$\mu_{\omega}$ | 3.5 | 3.472(0.028) | 3.475(0.028) | 3.469(0.028)
$\phi_{\theta,\theta}$ | 0.4 | 0.14(0.20) | 0.23(0.23) | 0.16(0.18)
$log\sigma_{\theta}$ | 2.996 | 2.77(0.15) | 2.88(0.19) | 2.70(0.14)
$\phi_{\theta,\theta}=0.4$, $tmax=25$, $n=100$
Parameters | True value | MLE | Prior1 | Prior2
no. divergent transitions | NA | NA | 2 | 0
$\mu_{\alpha}$ | 2 | 1.995(0.011) | 1.995(0.012) | 1.995(0.012)
$\mu_{\theta}$ | 20 | 20.2(8.7) | 18.3(17.5) | 20.1(7.4)
$\mu_{\omega}$ | 3.5 | 3.506(0.022) | 3.508(0.022) | 3.504(0.021)
$\phi_{\theta,\theta}$ | 0.4 | 0.50(0.17) | 0.59(0.18) | 0.46(0.13)
$log\sigma_{\theta}$ | 2.996 | 3.25(0.18) | 3.43(0.28) | 3.13(0.14)
$\phi_{\theta,\theta}=0.4$, $tmax=50$, $n=25$
Parameters | True value | MLE | Prior1 | Prior2
no. divergent transitions | NA | NA | 0 | 0
$\mu_{\alpha}$ | 2 | 1.974(0.015) | 1.974(0.015) | 1.973(0.015)
$\mu_{\theta}$ | 20 | 20.0(3.8) | 19.8(4.9) | 20.1(4.2)
$\mu_{\omega}$ | 3.5 | 3.520(0.032) | 3.523(0.032) | 3.515(0.031)
$\phi_{\theta,\theta}$ | 0.4 | 0.42(0.14) | 0.48(0.15) | 0.42(0.13)
$log\sigma_{\theta}$ | 2.996 | 2.84(0.13) | 2.92(0.16) | 2.79(0.13)
$\phi_{\theta,\theta}=0.4$, $tmax=50$, $n=100$
Parameters | True value | MLE | Prior1 | Prior2
no. divergent transitions | NA | NA | 0 | 0
$\mu_{\alpha}$ | 2 | 1.9865(0.0076) | 1.9864(0.0076) | 1.9861(0.0076)
$\mu_{\theta}$ | 20 | 20.7(3.9) | 20.0(5.0) | 20.7(4.1)
$\mu_{\omega}$ | 3.5 | 3.512(0.015) | 3.513(0.015) | 3.510(0.015)
$\phi_{\theta,\theta}$ | 0.4 | 0.41(0.13) | 0.47(0.15) | 0.41(0.12)
$log\sigma_{\theta}$ | 2.996 | 2.89(0.12) | 2.97(0.17) | 2.85(0.11)

From Table 1 we find that both frequentist and Bayesian inferences show good
estimates for $\mu_{\alpha}$ and $\mu_{\omega}$. It is interesting to see that
the auto-correlation for $\theta_{t}$ is not always under-estimated under all
settings (for example $(tmax=25,n=50)$), this can be also seen from the tables
for parameter estimates in Supporting Information. Bayesian inference with
Prior1 (half-Cauchy prior) generally reports smaller estimates of
$\mu_{\theta}$ than MLE and Prior2 (lognormal prior) but larger estimates of
$\phi_{\theta,\theta}$ and $log\sigma_{\theta}$. The estimates with MLE and
Prior2 are close to each other while the estimates with Prior2 show fewer
uncertainties for $\phi_{\theta,\theta}$ and $log\sigma_{\theta}$ implied by
the smaller standard errors in the brackets. Prior2 also reports smaller
estimates for $log\sigma_{\theta}$ compared with MLE and Prior1 since it puts
very large weight on small values of the variance, as will be graphically
demonstrated in section 3.4. We also find that $\phi_{\theta,\theta}$ and
$\log\sigma_{\theta}$ are difficult parameters to estimate since none of these
three techniques can estimate them accurately across all the cases. However,
the estimates are based on one realization of simulation, the discrepancy
between estimates to the true value would vary from simulation to simulation.

We also compare the estimates across the different sample sizes. We typically
compare the estimates between setting $(tmax=25,n=50)$ and $(tmax=25,n=100)$,
$(tmax=50,n=25)$ and $(tmax=50,n=100)$, $(tmax=25,n=100)$ and
$(tmax=50,n=100)$. We find that increasing the mean sample size at each time
point does not necessarily increase the certainty of the estimates, but the
data set with increased time points $(tmax=50,n=100)$ contains more
information on the parameters of interest and thus reports more certain
estimates compared with the data set with $(tmax=25,n=100)$. The same
conclusion can be also drawn by making similar comparisons among the estimates
in Table S1 and S2 in Supporting Information.

We can also find from Table 1, Table S1 and S2 from Supporting Information
that the Bayesian inference with Prior1 in some cases report 1 or 2 divergent
transitions while with Prior2 there are no divergent transitions reported.
This implies that the geometric shape of posterior likelihood with Prior1 is
more challenging for sampling probably due to light tails and thus potentially
leads to an incomplete exploration of the target distribution.

###  3.3 Bias Plot

Figure 1: Bias plots for the auto-regressive parameter $\phi_{\theta,\theta}$
(the upper plot) and for the scale parameter log$\sigma_{\theta}$ (the lower
plot) respectively under the setting with time series length $tmax=50$,
average annual sample size $n=25$, autocorrelation in $\theta_{t}$
$\phi_{\theta,\theta}=0.4$ and 44 replications (50 replications were
conducted, among which 6 replications report 3 or more divergent transitions
for the MCMC of Bayesian inference and thus are removed from the analysis).

The comparison between the estimates in the last section is based on one
realization of the simulation. To make comparisons of estimates over more
realizations, the simulation was repeated 50 times under the setting of
$(tmax=50,n=25)$. Due to divergent transitions, only 44 out of 50 replicates
were kept and the replications with more than 0.1% divergent transitions (in
2000 iterations) were excluded from the analysis. For the estimate of
$\phi_{\theta,\theta}$ and $log\sigma_{\theta}$ in each replication, the bias
was calculated in a frequentist framework as the absolute difference between
the true value and the mean estimate from each inference technique. The
absolute bias for $\phi_{\theta,\theta}$ and $log\sigma_{\theta}$ are
graphically displayed in the upper and lower plot in Fig. 1 respectively. From
the upper plot we find that in most replications, Bayesian inference with
Prior1 slightly outperforms the frequentist inference and Bayesian inference
with Prior2, the latter two reported very close estimates for
$\phi_{\theta,\theta}$. One striking thing is that the bias is close to or
even larger than 0.4 for some replications, this suggests that the inferences
report even negative estimates of $\phi_{\theta,\theta}$ and it again turns
out to be a difficult parameter. In the lower plot, we can see no single
inference technique stands out in estimating the scale parameter
$log\sigma_{\theta}$.

###  3.4 Prior-posterior distribution

(a) $(tmax=25,n=50)$.

(b) $(tmax=25,n=100)$.

(c) $(tmax=50,n=25)$.

(d) $(tmax=50,n=100)$.

Figure 2: Histograms of posterior samples of the scale parameter
$\sigma_{\theta}$ from models with two different prior distributions. Plot
(a), (b), (c) and (d) correspond to sample size setting $(tmax=25,n=50)$,
$(tmax=25,n=100)$, $(tmax=50,n=25)$ and $(tmax=50,n=100)$ respectively. On
each subplot, the left one shows the histogram of posterior samples given
half-Cauchy (0, 10) prior on $\sigma_{\theta}$ and similarly, the right one
displays the histogram of posterior samples given log-normal (1, 0.5) prior on
$\sigma_{\theta}$. Overlain on each subplot (the solid black lines) is the
corresponding prior density function. The red lines indicate the true value of
$\sigma_{\theta}$. Only $\phi_{\theta,\theta}=0.4$ was considered in the
simulations.

Fig. 2 shows histograms of posterior samples of the scale parameter
$\sigma_{\theta}$ from models with the two different prior distributions:
half-Cauchy (0, 10) and log-normal (1, 0.5), which are represented by solid
lines in the left and right plot on each subplot respectively. The true value
of $\sigma_{\theta}$ is indicated by a solid red line. Plot (a), (b), (c) and
(d) correspond to setting $(tmax=25,n=50)$, $(tmax=25,n=100)$,
$(tmax=50,n=25)$ and $(tmax=50,n=100)$ respectively. We can see from plot (a)
that the priors are quite informative and pull the posteriors towards small
values away from the true value and this prior-domination is more clear with
log-normal prior where the prior distribution sharply peaks at 2. The
domination is not mitigated even though the mean annual sample size is
increased to 100 as shown in plot (b). With the same total sample size in plot
(c) $(tmax=50,n=25)$ as that in plot (a) $(tmax=25,n=50)$, the posterior
likelihoods in plot (c) are, however, not dominated by the priors. The prior-
domination is also mitigated in plot (d) compared with plot (b) by increasing
the time points from 25 to 50.

Altogether, the informative log-normal prior pulls more of the posterior
towards a narrower range of smaller parameter values especially when the
number of time points in the data is small. The posterior samples are less
dominated by the half-Cauchy prior in this case. Increasing the annual mean
sample size does not necessarily lead to better identification of the small
region of parameter space. Only the amount of time points is the matter for
the likelihood to overwhelm the prior distribution and to dominate the
posterior distribution.

###  3.5 Computational efficiency with and without Laplace approximation

In tmbstan, sampling can be performed with or without Laplace approximation
for the random effects. It is possible to mix the Laplace approximation with
MCMC by specifying laplace=TRUE, such that the random effects are integrated
with the Laplace approximation in TMB and other parameters (such as fixed
effects and hyperparameters specifying the distribution of the random effects)
are handled by the NUTS in Stan. In the case studies in monnahan2018no, the
Bayesian inference algorithms with Laplace approximation is less
computationally efficient than without Laplace approximation, where the
efficiency is defined as the minimum effective sample size per second.
Following that definition, we calculated the efficiency of tmbstan with and
without Laplace approximation with simulated data. Different from
monnahan2018no, we did not consider the computational efficiency of
Frequentist inference with the Laplace approximation, as explained in the last
section.

In Fig. 3, plot (a) displays violin plots of computational efficiency without
(orange) and with (green) Laplace approximation (la) of Bayesian inference
with Prior1 under different sample size settings. The setting 1, 2, 3, 4 on x
axis stand for setting $(tmax=25,n=50)$, $(tmax=25,n=100)$, $(tmax=50,n=25)$
and $(tmax=50,n=100)$ respectively. Only $\phi_{\theta,\theta}=0.4$ was
considered and the divergent transitions were not taken into account. Inside
the violin plots are box plots showing the quantiles of 50 realized
computational efficiencies. Similarly, the violin plots of computational
efficiency with Prior2 are shown on plot (b). We find from both plot (a) and
(b) that Bayesian inference without Laplace approximation generally is more
efficient under setting 1, 2, and 3, the outperformance is more manifest when
the sample size is small $(tmax=25,n=50)$. However, when the sample size is
increased to $(tmax=50,n=100)$, inference with Laplace approximation turns out
to be slightly more efficient than that without Laplace approximation, the
boxplots and violin plots also tend to be more compact under this setting.

(a) Computational efficiency with Prior1.

(b) Computational efficiency with Prior2.

Figure 3: Violin plots of computational efficiency (minimum effective sample
size per second) without (orange) and with (green) Laplace approximation (la).
The four settings on x axis correspond to sample size setting
$(tmax=25,n=50)$, $(tmax=25,n=100)$, $(tmax=50,n=25)$ and $(tmax=50,n=100)$
respectively. Plot (a) shows the computational efficiency of Bayesian
inference with Prior1 and plot (b) with Prior2. Only
$\phi_{\theta,\theta}=0.4$ was used in simulations. Inside the violin plots
are box plots showing the quantiles of 50 realized computational efficiencies.
For each realization among the 50 simulations and across the settings, the
same specifications in tmbstan are used.

Even though the technique in which the random effects are integrated out by
Laplace approximation in TMB turns out to be less efficient in most settings,
we still provide a counterexample from monnahan2018no in which the enabling of
Laplace approximation is always less computationally efficient in the case
studies.

###  3.6 Laplace approximation check

(a) $(tmax=25,n=50)$ with Prior2.

(b) $(tmax=25,n=100)$ with Prior2.

(c) $(tmax=50,n=25)$ with Prior2.

(d) $(tmax=50,n=100)$ with Prior2.

Figure 4: Pair plots of posterior samples for Laplace approximation check from
one realization of the simulation with Prior2. The four plots (a) (b) (c) and
(d) correspond to the four settings of sample size in simulation. The random
effects in the TMB model can be integrated with two techniques: (1) full MCMC
integration via NUTS and (2) Laplace approximation. To check the accuracy of
Laplace approximation to the posterior likelihood density, the posterior
samples for all the fixed effects in the model without (yellow dots) and with
Laplace approximation (green dots) are shown pair-wisely on the same plot.
Columns and rows on the lower diagonal correspond to pair-wise parameters,
with the diagonal showing QQ-plot of posterior samples from Bayesian inference
without (yellow dots) and with (green dots) Laplace approximation for that
parameter including a 1:1 line in yellow. The large red circles on the off-
diagonal plots represent the pairwise means. On each off-diagonal plot, there
are 4000 yellow dots corresponding to 1000 samples retained from each of four
chains without Laplace approximation, so as the green dots with Laplace
approximation. Posterior rows were randomized to prevent consistent
overplotting of one integration technique. Overlaps in the two colored dots
suggest that the Laplace approximation is accurate.

By comparing the Bayesian posteriors with and without Laplace approximation,
we are allowed to check how well the Laplace approximation works. Fig. 4 shows
pair plots of posterior samples with and without Laplace approximation done by
TMB under different sample size settings with Prior2. Only autocorrelation in
$\theta_{t}$ $\phi_{\theta,\theta}=0.4$ was considered. Plot (a), (b), (c) and
(d) correspond to setting $(tmax=25,n=50)$, $(tmax=25,n=100)$,
$(tmax=50,n=25)$ and $(tmax=50,n=100)$ respectively. On each subplot, the
lower diagonal plots contain pairwise parameter posterior points. The green
dots represent posterior points from full MCMC integration via NUTS and the
yellow points from enabled Laplace approximation of the random effects. The
hollow red circles on the off-diagonal plots represent the pairwise means. The
diagonal shows QQ-plot of posterior samples from Bayesian inference without
(yellow dots) and with (green dots) Laplace approximation for that parameter
including a 1:1 line in yellow. Even though the posterior points are densely
packed, the overlap of the red circles with each technique shows seemingly
good alignment of the two versions of the posterior, and this suggests that
the Laplace approximation to the marginal likelihood where random effects are
integrated out works well. Similar pair plots for Laplace approximation check
with Prior1 can be found in Supporting Information.

##  4 Real-data case study

Having established the utility of our modeling approach and frequentist and
Bayesian inference in the context of simulated data, we also applied the same
statistical model to the analysis of a real great tit dataset of practical
interest. The observed data were collected from a Dutch great tit (Parus
major) population at the Hoge Veluwe National Park in the Netherlands (52°02’
- 52°07’N, 5°51’ - 5°32E). The recorded variables include the number of
chicks, number of fledglings, mother ID, brood laying date and so on for each
brood. Laying dates are presented as the number of days after March 31 (day
1=April 1, day 31=May 1). Similar to Reed2013, only the broods with one or
more chicks were considered in our analysis due to the high proportion (15.7%)
of zero-observations in the number of fledglings among the broods. The number
of fledglings was taken as the fitness component and assumed to be Poisson
distributed. The analyzed dataset consists of brood records breeding in 61
years from 1955 to 2015 and the sample size in a specific year ranges from 10
to 164 with an average of 81 across the study years. See Reed2013 for more
details on the study population and fieldwork procedures.

The focus of this empirical study is to compare the computational efficiency
of Bayesian inference with and without Laplace approximation and to check the
accuracy of Laplace approximation. However, since the true structure of the
model is unknown, we first conducted model selection under the frequentist
framework and the candidate models considered are different from each other
only in the model structure of stochastic $\alpha_{t}$, $\theta_{t}$ and
$\omega_{t}$. The details of all the candidate models including the best model
are given in Supporting Information. We then made Bayesian inference with the
two different priors as in the simulation study using the selected model. For
each prior distribution, we implemented tmbstan with and without Laplace
approximation to check the accuracy of Laplace approximation.

Table 2 lists the reported estimates of model parameters from maximum
likelihood (MLE) and Bayesian estimates with half-Cauchy (0, 10) prior
(Prior1) and log-normal (1, 0.5) prior (Prior2). The best model indicates
VAR(1) structure of $\alpha_{t}$ and $\theta_{t}$ and non-zero correlation
$\hat{\rho}_{\alpha,\theta}$. The width of stabilizing fitness function turned
to be constant over the study years implied by zero $\hat{\omega}_{t}$.
Frequentist inference and Bayesian inference with Prior2 report close
estimates for $\phi_{\theta,\theta}$ but the estimates with Prior2 show again
less uncertainty for most of the estimates except for $\rho_{\alpha,\theta}$.
In terms of $\log\sigma_{\theta}$, Bayesian inference with Prior1 reports the
largest estimate and least certainty compared with the other two techniques.
The close resemblance between estimates of $\log\sigma_{\theta}$ based on
maximum likelihood and Bayesian inferences suggests that the data contains a
good amount of information on $\log\sigma_{\theta}$ so that the maximum
likelihood overwhelms the log-normal prior and dominates the posterior
likelihood.

Table 2: Frequentist and Bayesian estimates of parameters in the selected model with great tit dataset. The Bayesian estimates (in column Prior1 and Prior2) are obtained without Lapalace approximation done by TMB.  parameter | MLE | Prior1 | Prior2
---|---|---|---
$\mu_{\alpha}$ | 2(0.0369) | 2(0.0491) | 2(0.0379)
$\mu_{\theta}$ | 18.5(5.35) | 18.8(7.12) | 19.4(5.09)
$\mu_{\omega}$ | 3.88(0.055) | 3.89(0.0563) | 3.86(0.0522)
$\phi_{\alpha,\alpha}$ | 0.379(0.12) | 0.458(0.13) | 0.398(0.124)
$\phi_{\theta,\theta}$ | 0.48(0.112) | 0.545(0.114) | 0.477(0.102)
$log\sigma_{\alpha}$ | -1.72(0.14) | -1.63(0.152) | -1.76(0.126)
$log\sigma_{\theta}$ | 3.07(0.137) | 3.16(0.155) | 2.98(0.125)
$\rho_{\alpha,\theta}$ | -0.728(0.0825) | -0.715(0.0895) | -0.661(0.0987)

Table 3 shows computational efficiencies of Bayesian inference without and
with Laplace approximation. It turns out that the computational efficiency
with Laplace approximation is approximately half of that without Laplace
approximation in both models with Prior1 and Prior2.

Table 3: Comparison of computational efficiency between Bayesian inference without (in the row "Full MCMC") and with Laplace approximation (in the row "Laplace approximation") for random effects for the great tit case study. Model | Inference | Time(s) | min.ESS | Efficiency(ESS/t)
---|---|---|---|---
Prior 1 | Full MCMC | 1542.215 | 186.7651 | 0.1211019
Laplace approximation | 15491.85 | 1004.643 | 0.06484975
Prior 2 | Full MCMC | 1266.096 | 291.0717 | 0.229897
Laplace approximation | 7815.218 | 1111.257 | 0.1421914

Similar to Fig. 4, Fig. 5 and Fig. 6 display pair plots of posterior samples
to check the accuracy of Laplace approximation with Prior1 and Prior2
respectively. Both the figures seemingly suggest a good mix of posterior
samples with and without Laplace approximation for all the parameters in the
selected model, indicating that the Laplace approximation assumption is met.

Figure 5: Pair plots of posterior samples for Laplace approximation test for
the great tit case study with Prior1. The random effects in the great tit TMB
model can be integrated with two techniques: (1) full MCMC integration via
NUTS and (2) Laplace approximation. To check the accuracy of Laplace
approximation to the posterior likelihood density, the posterior samples for
all the fixed effects in the model without (yellow dots) and with Laplace
approximation (green dots) are shown pair-wisely on the same plot. Columns and
rows on the lower diagonal correspond to pair-wise parameters, with the
diagonal showing QQ-plot of posterior samples from Bayesian inference without
(yellow dots) and with (green dots) Laplace approximation for that parameter
including a 1:1 line in yellow. The large red circles of the off-diagonal
plots represent the pairwise means. On each off-diagonal plot, there are 4000
yellow dots corresponding to 1000 samples retained from each of four chains
without Laplace approximation, so as the green dots with Laplace
approximation. Posterior rows were randomized to prevent consistent
overplotting of one integration technique. Overlaps in the two colored dots
suggest the Laplace approximation assumption is met. Figure 6: Pair plots of
posterior samples for Laplace approximation test for the great tit case study
with Prior2. The random effects in the great tit TMB model can be integrated
with two techniques: (1) full MCMC integration via NUTS and (2) Laplace
approximation. To check the accuracy of Laplace approximation to the posterior
likelihood density, the posterior samples for all the fixed effects in the
model without (yellow dots) and with Laplace approximation (green dots) are
shown pair-wisely on the same plot. Columns and rows on the lower diagonal
correspond to pair-wise parameters, with the diagonal showing QQ-plot of
posterior samples from Bayesian inference without (yellow dots) and with
(green dots) Laplace approximation for that parameter including a 1:1 line in
yellow. The large red circles of the off-diagonal plots represent the pairwise
means. On each off-diagonal plot, 4000 yellow dots correspond to 1000 samples
retained from each of four chains without Laplace approximation, so as the
green dots with Laplace approximation. Posterior rows were randomized to
prevent consistent overplotting of one integration technique. Overlaps in the
two colored dots suggest the Laplace approximation assumption is met.

##  5 Conclusions and extensions

In this study, we have investigated frequentist inference and Bayesian
inference with two different priors. The inferences were implemented with a
state-space model estimating temporal fluctuating selection and with simulated
biological data under four different simulation settings. A state-of-the-art R
package (tmbstan) for fast fitting statistical models was used for Bayesian
inference with Laplace approximation turning on or off. The simulation studies
show that the choice of prior can have an important impact on the geometric
shape of the posterior distributions of the model parameters and a non-
informative prior (in this study uniform prior and inverse-gamma prior on the
scale parameter) may lead to unstable inference since the Markov chains may
not converge or get stuck in part of the ridge of posterior. With unobserved
states following a VAR(1) process, we also found that the autoregressive
parameters and the scale parameters in the variance-covariance matrix of the
states are difficult and challenging to be estimated accurately. The increased
sample size at each time point does not necessarily provide more information
for the transition parameters and scale parameters. Only more time points in
the data could make the likelihood dominate the posterior likelihood and thus
lead to better estimates of these parameters. Half-Cauchy prior on the scale
parameter leads to less stable inference than log-normal prior indicated by
the number of divergent transitions in the Markov Chains. Laplace
approximation for the random effects turns out to be accurate suggested by the
pair plots of the posterior samples with and without Laplace approximation for
both the simulation studies and the great tit case study. Turning on Laplace
approximation in tmbstan would probably reduce computational efficiency but it
is worth trying when there is a good amount of data, in which case the Laplace
approximation is more likely to be accurate and also potentially improve the
computational efficiency of MCMC.

In our study, we used arbitrary prior distributions, however, the prior
information can be obtained from different sources. For example, in our great
tit case study, the timing and width of the caterpillar peak can provide a
clue for the time window of optimal laying dates, thus the information can be
used to decide the prior for the scale parameter of the optimal laying dates.
Prior information can also be generated from previous studies on the same
species and more general ecological knowledge coming from other related
species Tufto . (2000).

We conducted simulation studies with only AR(1) process of the optimal laying
dates, but the model is formulated and coded in a way that can be effortlessly
extended to order-1 vector autoregression (VAR(1)). It can be widely used for
modeling ecological processes where auto-correlation and cross-correlation in
the processes arise due to shared environmental variables at either temporal
or spatial scale. We expect more ecologists to adopt these two new estimation
methods, TMB, and tmbstan, given its flexibility in either frequentist or
Bayesian inference for a wide range of models, including the models where the
unobserved ecological processes are treated as latent variables and assumed to
be VAR processes. However, the drawback of Bayesian VAR (BVAR) methods is that
it usually requires estimation of a large number of parameters and thus the
over-parameterization might lead to unstable inference and inaccurate out-of-
sample forecasts. Some shrinkage methods (Sims Zha, 1998; Koop ., 2010;
Giannone ., 2015; Sørbye Rue, 2017, for example) were thereby developed, in
which Bayesian priors provide a logical and consistent method of imposing
parameter restrictions that can be potentially applied to ecological data
cases.

## References

  * Albertsen . (2016) albertsen2016choosingAlbertsen, CM., Nielsen, A. Thygesen, UH.  2016\.  Choosing the observational likelihood in state-space stock assessment models Choosing the observational likelihood in state-space stock assessment models.  Canadian Journal of Fisheries and Aquatic Sciences745779–789.  10.1139/cjfas-2015-0532
  * Auger-Méthé . (2017) auger2017spatiotemporalAuger-Méthé, M., Albertsen, CM., Jonsen, ID., Derocher, AE., Lidgard, DC., Studholme, KR.Flemming, JM.  2017\.  Spatiotemporal modelling of marine movement data using Template Model Builder (TMB) Spatiotemporal modelling of marine movement data using Template Model Builder (TMB).  Marine Ecology Progress Series565237–249.  10.3354/meps12019
  * Betancourt Girolami (2015) betancourt2015hamiltonianBetancourt, M. Girolami, M.  2015\.  Hamiltonian Monte Carlo for hierarchical models Hamiltonian Monte Carlo for hierarchical models.  Current trends in Bayesian methodology with applications7930.
  * Cadigan (2015) cadigan2015stateCadigan, NG.  2015\.  A state-space stock assessment model for northern cod, including under-reported catches and variable natural mortality rates A state-space stock assessment model for northern cod, including under-reported catches and variable natural mortality rates.  Canadian Journal of Fisheries and Aquatic Sciences732296–308.  10.1139/cjfas-2015-0047
  * Cao . (2019) cao2019timeCao, Y., Visser, ME. Tufto, J.  2019\.  A time-series model for estimating temporal variation in phenotypic selection on laying dates in a Dutch great tit population A time-series model for estimating temporal variation in phenotypic selection on laying dates in a dutch great tit population.  Methods in Ecology and Evolution1091401–1411.
  * Carpenter . (2017) carpenter2017stanCarpenter, B., Gelman, A., Hoffman, MD., Lee, D., Goodrich, B., Betancourt, M.Riddell, A.  2017\.  Stan: A probabilistic programming language Stan: A probabilistic programming language.  Journal of statistical software761.
  * Chevin . (2015) Chevin2015Chevin, LM., Visser, ME. Tufto, J.  2015\.  Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection.  Evolution6992319–2332.  10.1111/evo.12741
  * Fuglstad . (2019) fuglstad2019intuitiveFuglstad, GA., Hem, IG., Knight, A., Rue, H. Riebler, A.  2019\.  Intuitive principle-based priors for attributing variance in additive model structures Intuitive principle-based priors for attributing variance in additive model structures.  arXiv preprint arXiv:1902.00242.
  * Gelman . (2015) gelman2015stanGelman, A., Lee, D. Guo, J.  2015\.  Stan: A probabilistic programming language for Bayesian inference and optimization Stan: A probabilistic programming language for bayesian inference and optimization.  Journal of Educational and Behavioral Statistics405530–543.
  * Gelman . (2006) gelman2006priorGelman, A. .  2006\.  Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper) Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).  Bayesian analysis13515–534.
  * Giannone . (2015) giannone2015priorGiannone, D., Lenza, M. Primiceri, GE.  2015\.  Prior selection for vector autoregressions Prior selection for vector autoregressions.  Review of Economics and Statistics972436–451.
  * Grant Grant (2002) grant2002unpredictableGrant, PR. Grant, BR.  2002\.  Unpredictable evolution in a 30-year study of Darwin’s finches Unpredictable evolution in a 30-year study of Darwin’s finches.  science2965568707–711.  10.1126/science.1070315
  * Hoffman Gelman (2014) hoffman2014noHoffman, MD. Gelman, A.  2014\.  The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research1511593–1623.
  * Jeffreys Jeffreys (1961) jeffreys1961theoryJeffreys, H. Jeffreys, H.  1961\.  Theory of Probability (3rd edn). Theory of probability (3rd edn).  Oxford.
  * Karakani . (2016) karakani2016bayesianKarakani, HM., van Niekerk, J. van Staden, P.  2016\.  Bayesian Analysis of AR (1) model Bayesian analysis of AR (1) model.  arXiv preprint arXiv:1611.08747.
  * Koop . (2010) koop2010bayesianKoop, G., Korobilis, D. .  2010\.  Bayesian multivariate time series methods for empirical macroeconomics Bayesian multivariate time series methods for empirical macroeconomics.  Foundations and Trends® in Econometrics34267–358.
  * Kristensen . (2016) kristensen2016tmbKristensen, K., Nielsen, A., Berg, CW., Skaug, H. Bell, BM.  2016\.  TMB: Automatic Differentiation and Laplace Approximation TMB: Automatic differentiation and Laplace approximation.  Journal of Statistical Software7051–21.  10.18637/jss.v070.i05
  * Lei . (2011) lei2011bayesianLei, G., Boys, R., Gillespie, C., Greenall, A. Wilkinson, D.  2011\.  Bayesian inference for sparse VAR (1) models, with application to time course microarray data Bayesian inference for sparse VAR (1) models, with application to time course microarray data.  Journal of Biometrics and Biostatistics.
  * Lunn . (2000) lunn2000winbugsLunn, DJ., Thomas, A., Best, N. Spiegelhalter, D.  2000\.  WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility.  Statistics and computing104325–337.
  * Monnahan Kristensen (2018) monnahan2018noMonnahan, CC. Kristensen, K.  2018\.  No-U-turn sampling for fast Bayesian inference in ADMB and TMB: Introducing the adnuts and tmbstan R packages No-U-turn sampling for fast bayesian inference in ADMB and TMB: Introducing the adnuts and tmbstan R packages.  PloS one135e0197954.
  * Monnahan . (2017) monnahan2017fasterMonnahan, CC., Thorson, JT. Branch, TA.  2017\.  Faster estimation of Bayesian models in ecology using Hamiltonian Monte Carlo Faster estimation of bayesian models in ecology using Hamiltonian Monte Carlo.  Methods in Ecology and Evolution83339–348.
  * Neal . (2011) neal2011mcmcNeal, RM. .  2011\.  MCMC using Hamiltonian dynamics MCMC using Hamiltonian dynamics.  Handbook of markov chain monte carlo2112.
  * Reed . (2013) Reed2013Reed, TE., Jenouvrier, S. Visser, ME.  2013\.  Phenological mismatch strongly affects individual fitness but not population demography in a woodland passerine Phenological mismatch strongly affects individual fitness but not population demography in a woodland passerine.  Journal of Animal Ecology821131–144.  10.1111/j.1365-2656.2012.02020.x
  * Schlaifer Raiffa (1961) raiffa1961appliedSchlaifer, R. Raiffa, H.  1961\.  Applied statistical decision theory Applied statistical decision theory.  Wiley Cambridge.
  * Sims Zha (1998) sims1998bayesianSims, CA. Zha, T.  1998\.  Bayesian methods for dynamic multivariate models Bayesian methods for dynamic multivariate models.  International Economic Review949–968.
  * Sørbye Rue (2017) sorbye2017penalisedSørbye, SH. Rue, H.  2017\.  Penalised complexity priors for stationary autoregressive processes Penalised complexity priors for stationary autoregressive processes.  Journal of Time Series Analysis386923–935.
  * Stan Development Team (20181) team2018rstanStan Development Team.  20181\.  RStan: the R interface to Stan Rstan: the R interface to stan.  R package version 2.17.3.  http://mc-stan.org.
  * Stan Development Team (20182) stan2018stanStan Development Team.  20182\.  Stan Modeling Language Users Guide and Reference Manual Stan modeling language users guide and reference manual.  Version 2.18.0.  http://mc-stan.org.
  * Tufto . (2012) tufto2012estimatingTufto, J., Lande, R., Ringsby, TH., Engen, S., Sæther, BE., Walla, TR. DeVries, PJ.  2012\.  Estimating Brownian motion dispersal rate, longevity and population density from spatially explicit mark–recapture data on tropical butterflies Estimating Brownian motion dispersal rate, longevity and population density from spatially explicit mark–recapture data on tropical butterflies.  Journal of Animal Ecology814756–769.
  * Tufto . (2000) tufto2000bayesianTufto, J., Sæther, BE., Engen, S., Arcese, P., Jerstad, K., Røstad, OW. Smith, JN.  2000\.  Bayesian meta-analysis of demographic parameters in three small, temperate passerines Bayesian meta-analysis of demographic parameters in three small, temperate passerines.  Oikos882273–281.
  * Vedder . (2013) vedder2013quantitativeVedder, O., Bouwhuis, S. Sheldon, BC.  2013\.  Quantitative assessment of the importance of phenotypic plasticity in adaptation to climate change in wild bird populations Quantitative assessment of the importance of phenotypic plasticity in adaptation to climate change in wild bird populations.  PLoS Biology117e1001605.  10.1371/journal.pbio.1001605
  * Wei (2006) wei2006timeWei, W.  2006\.  Time Series Analysis: Univariate and Multivariate Methods Time series analysis: Univariate and multivariate methods (2nd ).  Pearson Addison Wesley.
  * Zellner (1986) zellner1986assessingZellner, A.  1986\.  On assessing prior distributions and Bayesian regression analysis with g-prior distributions On assessing prior distributions and Bayesian regression analysis with g-prior distributions.  Bayesian inference and decision techniques.