scripts/test_stan_code.R

In lorecatta/DENVclimate:

devtools::load_all()

foi_covariates <- readRDS(file.path("output", "foi_data_cov_rescaled.rds"))

m_values <- readRDS(file.path("output", "central_estimates_from_Mordecai", "m_mean_DayTemp_const_term.rds"))

T_inf <- 6 # from Zika paper(T_inf = 1/human recovry rate)

R0_values <- foi_covariates$R0_1

R0_values_n <- length(R0_values)

temp_values <- foi_covariates$DayTemp_const_term

N <- 100

out <- matrix(0, nrow = R0_values_n, ncol = N)

for (i in 1:N) {

  pr_vals <- prior_values(1)

  c_a <- pr_vals$c_a
  T0_a <- pr_vals$T0_a
  Tm_a <- pr_vals$Tm_a
  c_b <- pr_vals$c_b
  T0_b <- pr_vals$T0_b
  Tm_b <- pr_vals$Tm_b
  c_c <- pr_vals$c_c
  T0_c <- pr_vals$T0_c
  Tm_c <- pr_vals$Tm_c
  c_lf <- pr_vals$c_lf
  T0_lf <- pr_vals$T0_lf
  Tm_lf <- pr_vals$Tm_lf
  c_PDR <- pr_vals$c_PDR
  T0_PDR <- pr_vals$T0_PDR
  Tm_PDR <- pr_vals$Tm_PDR

  message("c_a = ", c_a)
  message("T0_a = ", T0_a)
  message("Tm_a = ", Tm_a)
  message("c_b = ", c_b)
  message("T0_b = ", T0_b)
  message("Tm_b = ", Tm_b)
  message("c_c = ", c_c)
  message("T0_c = ", T0_c)
  message("Tm_c = ", Tm_c)
  message("c_lf = ", c_lf)
  message("T0_lf = ", T0_lf)
  message("Tm_lf = ", Tm_lf)
  message("c_PDR = ", c_PDR)
  message("T0_PDR = ", T0_PDR)
  message("Tm_PDR = ", Tm_PDR)

  out[,i] <- R_0(temp_values,
                 m_values,
                 T_inf,
                 c_a, T0_a, Tm_a, c_b, T0_b, Tm_b, c_c, T0_c,
                 Tm_c, c_lf, T0_lf, Tm_lf, c_PDR, T0_PDR, Tm_PDR)

}