DENVclimate:

devtools::load_all()

library(reshape2)
library(ggplot2)


# define parameters -----------------------------------------------------------


parameters <- list(T_inf = 6, # from Zika paper(T_inf = 1/human recovry rate)
                   min_value = 0,
                   min_value_lf = 1,
                   ec = 1)

in_dir <- file.path("output", "central_estimates_from_Mordecai")

out_dir <- file.path("output", "central_estimates_from_Mordecai")


# load data -------------------------------------------------------------------


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

a_mean <- readRDS(file.path(in_dir, "a_mean_DayTemp_const_term.rds"))
b_mean <- readRDS(file.path(in_dir, "b_mean_DayTemp_const_term.rds"))
c_mean <- readRDS(file.path(in_dir, "c_mean_DayTemp_const_term.rds"))
PDR_mean <- readRDS(file.path(in_dir, "PDR_mean_DayTemp_const_term.rds"))
lf_mean <- readRDS(file.path(in_dir, "lf_mean_DayTemp_const_term.rds"))


# define variables ------------------------------------------------------------


R0_values <- foi_covariates$R0_1

temp_values <- foi_covariates$DayTemp_const_term


# extract traits estimated at observed T --------------------------------------


k_values <- a_mean[, "mean"]
beta_mh_values <- b_mean[, "mean"]
beta_hm_values <- c_mean[, "mean"]
PDR_values <- PDR_mean[, "mean"]
lf_values <- lf_mean[, "mean"]


# plot ------------------------------------------------------------------------


dat <- data.frame(temp = temp_values,
                  k = k_values,
                  beta_mh = beta_mh_values,
                  beta_hm = beta_hm_values,
                  PDR = PDR_values,
                  lf = lf_values)

dat_m <- melt(dat,
              id.vars = "temp",
              variable.name = "trait")

p <- ggplot(dat_m, aes(temp, value)) +
  geom_point() +
  facet_wrap(~ trait, scales = "free_y") +
  scale_x_continuous(expression("Temperature " (degree*C)),
                     limits = c(20, 50),
                     breaks = seq(20, 50, 5)) +
  scale_y_continuous("Trait") +
  theme_bw()

ggsave(filename = file.path("figures", "traits_at_observed_T.png"), plot = p,
       width = 18, height = 11, units = "cm", dpi = 300)


# derive some traits ----------------------------------------------------------

min_value <- parameters$min_value
min_value_lf <- parameters$min_value_lf

k_values[k_values < min_value] <- min_value
beta_mh_values[beta_mh_values < min_value] <- min_value
beta_hm_values[beta_hm_values < min_value] <- min_value
lf_values[lf_values < min_value_lf] <- min_value_lf

mu_values <- 1 / lf_values
eip_values <- 1 / PDR_mean[, "mean"]

m_values <- calculate_m(R0_values,
                        k_values,
                        beta_mh_values,
                        beta_hm_values,
                        mu_values,
                        eip_values,
                        parameters)


# plot ------------------------------------------------------------------------


dat_2 <- data.frame(temp = temp_values,
                    k = k_values,
                    beta_mh = beta_mh_values,
                    beta_hm = beta_hm_values,
                    PDR = PDR_values,
                    eip = eip_values,
                    lf = lf_values,
                    mu = mu_values,
                    m = m_values)

dat_m_2 <- melt(dat_2,
                id.vars = "temp",
                variable.name = "trait")

p <- ggplot(dat_m_2, aes(temp, value)) +
  geom_point() +
  facet_wrap(~ trait, scales = "free_y") +
  scale_x_continuous(expression("Temperature " (degree*C)),
                     limits = c(20, 50),
                     breaks = seq(20, 50, 5)) +
  scale_y_continuous("Trait") +
  theme_bw()

ggsave(filename = file.path("figures", "traits_at_observed_T_derived.png"), plot = p,
       width = 18, height = 15, units = "cm", dpi = 300)


# save m values ---------------------------------------------------------------


write_out_rds(m_values, out_dir, "m_mean_DayTemp_const_term")