DENVclimate:

devtools::load_all()

library(dplyr)
library(reshape2)
library(ggplot2)


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


in_dir <- file.path("output", "termal_response_fits", "informative")

params <- c("T0", "Tm", "c")


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


a_samps <- readRDS(file.path(in_dir, "a_samps.rds"))
b_samps <- readRDS(file.path(in_dir, "b_samps.rds"))
c_samps <- readRDS(file.path(in_dir, "c_samps.rds"))
PDR_samps <- readRDS(file.path(in_dir, "PDR_samps.rds"))
lf_samps <- readRDS(file.path(in_dir, "lf_DENV_samps.rds"))


# pre processing --------------------------------------------------------------


ids <- seq_len(nrow(a_samps))

a_samps_2 <- a_samps %>%
  mutate(trait = "a",
         id = ids)

b_samps_2 <- b_samps %>%
  mutate(trait = "b",
         id = ids)

c_samps_2 <- c_samps %>%
  mutate(trait = "c",
         id = ids)

PDR_samps_2 <- PDR_samps %>%
  mutate(trait = "PDR",
         id = ids)

lf_samps_2 <- lf_samps %>%
  mutate(trait = "lf",
         id = ids) %>%
  rename(c = qd) %>%
  select(-n.qd)

all_samps <- rbind(a_samps_2, b_samps_2, c_samps_2, PDR_samps_2, lf_samps_2)

all_samps_m <- melt(all_samps,
                    id.vars = c("id", "trait"),
                    measure.vars = params,
                    variable.name = "parameter")

all_samps_m$trait_par <- paste(all_samps_m$parameter, all_samps_m$trait, sep = "_")

all_samps_m$trait_par <- factor(all_samps_m$trait_par,
                                levels = c("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"))


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


p <- ggplot(data = all_samps_m, mapping = aes(x = value)) +
  geom_histogram(col = "white") +
  facet_wrap(~ trait_par, scales = "free") +
  scale_x_continuous("Trait value", labels = function(x) format(x, scientific = FALSE)) +
  scale_y_continuous("Frequency") +
  ggtitle("Informative priors (Posterior densities from Mordecai et al. 2017)") +
  theme_bw() +
  theme(axis.text = element_text(size = 6))

save_plot(p, "figures", "informative_priors", wdt = 20, hgt = 18)