## ----global_options, include=FALSE-----------------------------------------------------------------
knitr::opts_chunk$set(fig.align = 'center', warning=FALSE, message=FALSE, fig.asp=0.625, dev='png', global.par = TRUE, dev.args=list(pointsize=10), fig.path = 'figs/',  eval=TRUE, cache=TRUE)
library(MASS)

## ----setup, include=FALSE--------------------------------------------------------------------------
library(knitr)
local({
    hook_plot = knit_hooks$get('plot')
    knit_hooks$set(plot = function(x, options) {
        paste0('\n\n----\n\n', hook_plot(x, options))
    })
})


## ----eval=TRUE,echo=TRUE, warning=FALSE, results='hide',message=FALSE------------------------------
library(rstan)
library(dplyr)
library(lubridate)
library(ggplot2)
library(bayesplot)

theme_set(bayesplot::theme_default())
rstan_options(auto_write=TRUE)


## --------------------------------------------------------------------------------------------------
stan_dat_hier <- readRDS('pest_data_longer_stan_dat.RDS')


## ----comp-NB-hier-slopes, cache=TRUE, results="hide", message=FALSE--------------------------------
comp_model_NB_hier_slopes <- stan_model('hier_NB_regression_ncp_slopes_mod.stan')


## ----run-NB-hier-slopes----------------------------------------------------------------------------
fitted_model_NB_hier_slopes <-
  sampling(
    comp_model_NB_hier_slopes,
    data = stan_dat_hier,
    chains = 4#, #cores = 4,
#    control = list(adapt_delta = 0.95)
  )


## --------------------------------------------------------------------------------------------------
mcmc_hist(
  as.matrix(fitted_model_NB_hier_slopes, pars = "sigma_kappa"),
  binwidth = 0.005
)


## --------------------------------------------------------------------------------------------------
print(fitted_model_NB_hier_slopes, pars = c('kappa','beta','alpha','phi','sigma_mu','sigma_kappa','mu'))


## --------------------------------------------------------------------------------------------------
mcmc_hist(
  as.matrix(fitted_model_NB_hier_slopes, pars = "beta"),
  binwidth = 0.005
)


## ----ppc-full-hier-slopes--------------------------------------------------------------------------
y_rep <- as.matrix(fitted_model_NB_hier_slopes, pars = "y_rep")
ppc_dens_overlay(
  y = stan_dat_hier$complaints,
  yrep = y_rep[1:200,]
)


## ----ppc-group_max-hier-slopes-mean-by-mo----------------------------------------------------------
select_vec <- which(stan_dat_hier$mo_idx %in% 1:12)
ppc_stat_grouped(
  y = stan_dat_hier$complaints[select_vec],
  yrep = y_rep[,select_vec],
  group = stan_dat_hier$mo_idx[select_vec],
  stat = 'mean'
) + xlim(0, 11)


## ----comp-NB-hier-mos, cache=TRUE, results="hide", message=FALSE-----------------------------------
comp_model_NB_hier_mos <- stan_model('hier_NB_regression_ncp_slopes_mod_mos.stan')


## ----run-NB-hier-slopes-mos------------------------------------------------------------------------
fitted_model_NB_hier_mos <- sampling(comp_model_NB_hier_mos, data = stan_dat_hier, chains = 4, 
                                     #cores = 4, 
                                     control = list(adapt_delta = 0.9))


## ----ppc-full-hier-mos-----------------------------------------------------------------------------
y_rep <- as.matrix(fitted_model_NB_hier_mos, pars = "y_rep")
ppc_dens_overlay(
  y = stan_dat_hier$complaints,
  yrep = y_rep[1:200,]
)


## --------------------------------------------------------------------------------------------------
select_vec <- which(stan_dat_hier$mo_idx %in% 1:12)
ppc_stat_grouped(
  y = stan_dat_hier$complaints[select_vec],
  yrep = y_rep[,select_vec],
  group = stan_dat_hier$mo_idx[select_vec],
  stat = 'mean'
)


## --------------------------------------------------------------------------------------------------
# overlay prior density curve on posterior draws
gen_rho_prior <- function(x) {
  alpha <- 10; beta <- 5
  a <- -1; c <- 1
  lp <- (alpha - 1) * log(x - a) +
        (beta - 1) * log(c - x) -
        (alpha + beta - 1) * log(c - a) -
         lbeta(alpha, beta)
  return(exp(lp))
}
mcmc_hist(as.matrix(fitted_model_NB_hier_mos, pars = "rho"),
          freq = FALSE, binwidth = 0.01) +
  overlay_function(fun = gen_rho_prior) +
  xlim(-1,1)

## --------------------------------------------------------------------------------------------------
print(fitted_model_NB_hier_mos, pars = c('rho','sigma_mu','sigma_kappa', 'mu', 'kappa'))


## --------------------------------------------------------------------------------------------------
ppc_intervals(
  y = stan_dat_hier$complaints,
  yrep = y_rep,
  x = stan_dat_hier$traps
) +
  labs(x = "Number of traps", y = "Number of complaints")


## --------------------------------------------------------------------------------------------------
comp_model_NB_hier_slopes_llik <- stan_model('hier_NB_regression_ncp_slopes_mod_loglik.stan')

fitted_model_NB_hier_slopes_llik <-
  sampling(
    comp_model_NB_hier_slopes_llik,
    data = stan_dat_hier,
    chains = 4, #cores = 4,
    control = list(adapt_delta = 0.95)
  )


## --------------------------------------------------------------------------------------------------
comp_model_NB_hier_mos_llik <- stan_model('hier_NB_regression_ncp_slopes_mod_mos_loglik.stan')

fitted_model_NB_hier_mos_llik <-
  sampling(
    comp_model_NB_hier_mos_llik,
    data = stan_dat_hier,
    chains = 4, #cores = 4,
    control = list(adapt_delta = 0.95)
  )




## ----run-NB-hier-slopes-mos-loo--------------------------------------------------------------------
library(loo)
log_lik_slopes <- extract_log_lik(fitted_model_NB_hier_slopes_llik)
loo_slopes <- loo(log_lik_slopes, 
                  r_eff=relative_eff(exp(log_lik_slopes),rep(c(1,2,3,4),each=1000)))
loo_slopes
log_lik_mos <- extract_log_lik(fitted_model_NB_hier_mos_llik)
loo_mos <- loo(log_lik_mos,
               r_eff=relative_eff(exp(log_lik_mos),rep(c(1,2,3,4),each=1000)))
loo_mos

loo_compare(loo_slopes, loo_mos)