## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(dplyr)
library(ggplot2)
library(latex2exp)
library(bennu)

## ----run experiments, eval=FALSE----------------------------------------------
#  experiments <- expand_grid(
#    sigma = seq(0.05, 0.5, by = 0.05),
#    zeta = seq(0.05, 0.5, by = 0.05)
#  )
#  
#  experimental_validation_data <- tibble()
#  
#  for (i in 1:nrow(experiments)) {
#    sigma <- as.numeric(experiments[i, "sigma"])
#    zeta <- as.numeric(experiments[i, "zeta"])
#  
#    d <- model_random_walk_data(
#      zeta = zeta, sigma = sigma,
#      region_coeffs = c(5, 0.5, 1, 2),
#      c_region = c(-1, 2, 0, 1), N_t = 48
#    )
#    fit <- est_naloxone(d)
#  
#    row_results <- fit %>%
#      tidybayes::gather_draws(sigma, zeta) %>%
#      group_by(.variable) %>%
#      dplyr::summarise(
#        p50 = stats::quantile(.value, 0.5),
#        p25 = stats::quantile(.value, 0.25),
#        p75 = stats::quantile(.value, 0.75),
#        p05 = stats::quantile(.value, 0.05),
#        p95 = stats::quantile(.value, 0.95)
#      ) %>%
#      left_join(
#        tibble(
#          .variable = c("sigma", "zeta"),
#          true_value = c(sigma, zeta)
#        )
#      ) %>%
#      mutate(experiment = i)
#  
#    experimental_validation_data <- experimental_validation_data %>%
#      bind_rows(row_results)
#  }

## -----------------------------------------------------------------------------

p_res_sigma <- experimental_validation_data %>%
  group_by(experiment) %>%
  mutate(zeta_val = sum(true_value * as.numeric(.variable == "zeta"))) %>%
  ungroup() %>%
  filter(.variable == "sigma") %>%
  ggplot(aes(x=true_value,y=true_value)) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray20") +
  geom_errorbar(aes(y = p50,ymin=p05,ymax=p95, color = zeta_val,
                    group = experiment)) +
  geom_point(aes(y = p50, color = zeta_val)) +
  theme_classic() +
  labs(color = parse(text = TeX("$\\zeta$")),
       x = "True value", y = "Inferred value")


show(p_res_sigma)

## -----------------------------------------------------------------------------
p_res_zeta <- experimental_validation_data %>%
  group_by(experiment) %>%
  mutate(sigma_val = sum(true_value * as.numeric(.variable == "sigma"))) %>%
  ungroup() %>%
  filter(.variable == "zeta") %>%
  ggplot(aes(x=true_value,y=true_value)) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray20") +
  geom_errorbar(aes(y = p50,ymin=p05,ymax=p95, color = sigma_val,
                    group = experiment)) +
  geom_point(aes(y = p50, color = sigma_val)) +
  theme_classic() +
  labs(color = parse(text = TeX("$\\sigma$")),
       x = "True value", y = "Inferred value")

show(p_res_zeta)

## ----missing data experiments, eval=FALSE-------------------------------------
#  set.seed(43)
#  
#  rstan::rstan_options(auto_write = TRUE)
#  options(mc.cores = 4)
#  
#  
#  experiments <- expand_grid(reporting_freq = seq(1, 10, by = 1))
#  
#  missing_data_validation <- tibble()
#  
#  
#  for (ind in 1:nrow(experiments)) {
#    reporting_freq <- as.numeric(experiments[ind, "reporting_freq"])
#  
#    # generate complete data
#    d_complete <- bennu::model_random_walk_data(
#      zeta = 0.25, sigma = 0.25,
#      region_coeffs = c(5, 0.5, 1, 2),
#      c_region = c(-1, 2, 0, 1), N_t = 48
#    )
#  
#    # generate partial data based on reporting frequency
#    d_reported <- d_complete %>%
#      mutate(
#        nonreporting_times = times %% reporting_freq != 1,
#        Reported_Distributed = if_else(
#          times %% reporting_freq != 1,
#          NA_real_, Reported_Distributed
#        ),
#        Reported_Used = if_else(
#          times %% reporting_freq != 1,
#          NA_real_, Reported_Used
#        )
#      )
#  
#    if (reporting_freq == 1) {
#      d_reported <- d_complete
#    }
#  
#    fit <- est_naloxone(d_reported)
#  
#    row_used <- fit %>%
#      tidybayes::spread_draws(sim_used[i]) %>%
#      dplyr::left_join(
#        dplyr::mutate(d_complete, i = dplyr::row_number()),
#        by = "i"
#      ) %>%
#      group_by(.chain, .iteration, .draw) %>%
#      dplyr::summarise(
#        Actual_Used = sum(Orders * p_use),
#        Sim_Used = sum(sim_used)
#      ) %>%
#      ungroup() %>%
#      dplyr::mutate(p_diff = (Actual_Used - Sim_Used) / Actual_Used) %>%
#      dplyr::summarise(
#        p50 = stats::quantile(p_diff, 0.5),
#        p25 = stats::quantile(p_diff, 0.25),
#        p75 = stats::quantile(p_diff, 0.75),
#        p05 = stats::quantile(p_diff, 0.05),
#        p95 = stats::quantile(p_diff, 0.95)
#      ) %>%
#      mutate(reporting_freq = reporting_freq)
#  
#    missing_data_validation <- missing_data_validation %>% bind_rows(row_used)
#  }
#  

## ----missing data frequency plot----------------------------------------------

p_res_freq <- missing_data_validation %>%
  ggplot(aes(x=reporting_freq)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_errorbar(aes(y = p50,ymin=p05,ymax=p95), color = "#4876FF") +
  geom_errorbar(aes(y = p50,ymin=p25,ymax=p75), color = "#436EEE") +
  geom_point(aes(y = p50), color = "#27408B") +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "Frequency of reporting (months)", y = "Relative error")


show(p_res_freq)