Bayesian Estimation
All the parameters can be estimated either under a Poisson model with
argument family = "poisson" or under a Negative-Binomial
model which includes an additional overdispersion parameter \(\phi\) with argument
family="nb": \[
D_{xt} \sim Poisson (e_{xt}\mu_{xt}) \text{ or } D_{xt} \sim NB
(e_{xt}\mu_{xt},\phi)
\] Given the matrix of deaths deathFR and the matrix
of central exposures exposureFR, we can infer the posterior
distribution of all the parameters and obtain death rates forecasts for
the next 10 years under a Poisson model by a simple call to the
lc_stan function:
fitLC=lc_stan(death = deathFR,exposure=exposureFR, forecast = 10, family = "poisson",chains=1,iter=1000,cores=1)##
## SAMPLING FOR MODEL 'leecarter' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000235 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.35 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 3.77 seconds (Warm-up)
## Chain 1: 2.354 seconds (Sampling)
## Chain 1: 6.124 seconds (Total)
## Chain 1:Internally, the sampling is performed through the
rstan::sampling function. By default, Stan samples four
Markov chains of 2000 iterations. For each chain, there are 1000 warmup
iterations (hence, 4000 post warm-up draws in total). Moreover, by
default, Stan uses 1 core but we recommend using as many processors as
the hardware and RAM allow (up to the number of chains). All this can be
set using the following additional arguments : for instance, if one
wants to sample 2 chains of 3000 iterations with 1000 warm-up samples,
one can add chains=2,iter=3000,warmup=1000.
For additional arguments, the reader may refer to the
rstan::sampling documentation.
The output is an object of class stanfit (cf. the
rstan package) which contains, among others,
posterior draws, posterior summary statistics and convergence
diagnostics.
The easiest way to extract the posterior draws is to call the
rstan::extract function which returns a list with named
components corresponding to the model parameters.
## [1] "a" "b" "c" "ks" "sigma" "k" "phi"
## [8] "k_p" "mufor" "log_lik" "pos" "pos2" "pos3" "lp__"Among these parameters, we find
a: \(\alpha_x\).b: \(\beta_x\).k: \(\kappa_t\).k_p: forecasts of \(\kappa_t\).mufor: forecasts of death rates \(\mu_{xt}\).log_lik: log-likelihoods.phi: overdispersion parameter for the NB model.candsigma: drift and standard deviation of the random walk with drift.
The user can have access to an interface for interactive MCMC
diagnostics, plots and tables helpful for analyzing posterior samples
through the shinystan package (fore more details, you can
check the shinystan web
page)
Note: you need to close the shiny app before continuing in R.
Moreover, the package also includes functions to represent boxplots
of the posterior distribution based on the ggplot2
package.
For instance, boxplots for the parameters from the Lee-Carter model can be obtained via:
