Chapter 3 Simulation study: Proportion travelers

3.1 Bootstraping

f <- SEROPOSITIVE ~ work_out*edad + work_out*nm_sex + edu_cat + hist_fever + comm
R <- 999 # N bootstraps - 5 for testing, chage upon production

library(boot)
library(purrr)
library(furrr) #UPDATE
future::plan(multiprocess) #UPDATE

dat_sim <- dat %>%
  nest() %>%
  mutate(area = "Overall") %>%
  bind_rows(dat %>%
              mutate(area = ifelse(area == "0_periurban",
                                   "Periurban",
                                   "Rural")) %>%
              group_by(area) %>%
              nest()) %>%
  crossing(p_out = seq(0, 1, by = .1)) %>%
  mutate(data = map(.x = data, 
                    .f = ~mutate(.x, comm = as.factor(as.character(comm)))),
         boot = future_map2(.x = data, .y =  p_out, #UPDATE
                     .f = ~boot(data = .x,
                                statistic = standardization_sim, 
                                treatment = work_out, 
                                outcome = SEROPOSITIVE, 
                                p_out = .y,
                                formula1 = f,
                                R = R
                                )),
         boot_tab = map(.x = boot,
                        .f = ~tab_sim(.x)))

saveRDS(dat_sim, "./_out/_rds/result_04-travelers.rds")

dat_sim_res <- dat_sim %>%
  select(area, p_out, boot_tab) %>%
  unnest()

3.2 Plots

OBS = observed, NE = No Exposed, FE = Exposed, SE = Simulated Exposure

dat_sim_res %>%
  filter(type != "OBS") %>%
  filter(type != "NE") %>%
  filter(type != "FE") %>%
  filter(type != "FE - NE") %>%
  filter(type != "FE - SE") %>%
  #filter(type != "SE - OBS") %>% #UPDATE
  mutate(type = case_when( #UPDATE
    type == "SE - NE" ~ "RD [SE - NE]",
    type == "SE - OBS" ~ "RD [SE - NC]",
    type == "SE" ~ "SE"),
    type = fct_rev(type)) %>% #UPDATE
  #filter(type != "Simulated Exposure - No Exposed") %>%
  ggplot(aes(x = p_out, y = mean, col = area)) +
  # geom_vline(data = dat_sim_res %>%
  #                  filter(p_out == 0 &
  #                           type == "OBS") %>%
  #              select(-type) %>%
  #              crossing(type = c("SE",
  #                                "SE - OBS",
  #                                "SE - NE")),
  #            aes(xintercept = mean),
  #            linetype = "dashed") +
  geom_hline(data = . %>%
               filter(type != "SE"),
             aes(yintercept = 0), col = "grey50") +
  geom_rect(data = . %>%
               filter(type != "SE"),
            aes(ymin = -Inf, ymax = Inf, xmin = -Inf, xmax = Inf), #UPDATE
            alpha = .025, col = NA) +
  geom_linerange(aes(ymin = ll, ymax = ul), alpha = .6, 
                 size = 1.5) +
  geom_point(size = .7) +
  guides(col = F) +
  scale_color_npg() +
  #scale_color_viridis_c(option = "B", direction = -1) +
  facet_grid(type~area, scales = "free") +
  #geom_smooth(method = "lm", se = F, size = .2) +
  labs(y = "Standardized mean outcome", x = "Simulated proportion of inhabitants with outside village activities ") +
  theme_bw() +
  theme(legend.position = "top",
        #panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

ggsave("./_out/fig4.png", width = 7, height = 7,
       dpi = "retina", bg = "white")