In cmhoove14/DDNTD: Density Dependence in Neglected Tropical Diseases

knitr::opts_chunk$set(echo = TRUE)

require(tidyverse)

devtools::load_all()

sen_dat <- readRDS("../data/senegal_comms_village_level_summaries.rds")

ctrl_vils <- c("DF", "DT", "MB", "MT", "NM", "MG", "NE", "ST") #Communities in control group that only received MDA

Community-level summaries of Senegal communities

sen_dat %>% 
  filter(var %in% c("w", "post")) %>% 
  mutate(date = case_when(var == "w" ~ as.Date(paste0(year, "-03-01")),
                          var == "post" ~ as.Date(paste0(year, "-03-02")) + 1)) %>% 
  ggplot(aes(x = date, y = value, col = school)) +
    geom_line(size = 1.2) +
    theme_classic() +
    labs(x = "Time",
         y = "Mean worm burden estimate, (W)",
         title = "Worm burden over two years of MDA")

sen_dat %>% 
  filter(var == "Reff") %>% 
  ggplot(aes(x = school, y = value, fill = year, col = school)) +
    geom_bar(stat = "identity", 
             position = position_dodge(),
             width = 0.5,
             size = 1.2) +
    theme_classic() +
    theme(axis.text = element_text(size = 10),
          axis.title = element_text(size = 14),
          legend.position = "bottom") +
    scale_fill_manual(values = c("grey90", "grey10")) +
    ylim(c(0,10)) +
    labs(x = "Community",
         y = expression(italic(R[eff])~estimate),
         col = "Community")

sen_alphas <- sen_dat %>% 
  filter(var == "lambda" & year == 2017) %>% 
  mutate(alpha = value/(base_pars["omega"]*base_pars["theta"]*10))

comm_reff_profile <- function(comm){

  W_base = sen_dat %>% filter(school == comm & year == 2016 & var == "w") %>% pull(value)

  use_pars <- base_pars
  use_pars["alpha"] <- sen_alphas %>% filter(school == comm) %>% pull(alpha)
  fit_beta_Neq <- fit_pars_from_eq_vals(beta_guess = 2e-4, alpha_guess = 2e-4, Neq_guess = 1000, 
                                        W = W_base, Ip = 0.025, pars = use_pars)
  use_pars["beta"] <- fit_beta_Neq[1]

  reff_crv <- data.frame(W = exp_seq(1e-4, 100, 200)) %>% 
    mutate(Reff = sapply(W, Reff_W_kap, kap = 0.05, 
                         pars = use_pars)[1,],
           school = comm)

  return(reff_crv)

}

comm_reff_profiles <- bind_rows(lapply(unique(sen_dat$school), comm_reff_profile))

comm_reff_profiles %>% 
  ggplot(aes(x = W, y = Reff, col = school)) +
    geom_line(size = 1.2) +
    theme_classic() +
    scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-2,1,100)) +
    labs(x = "mean worm burden (W)",
         y = expression(italic(R[eff])))

comm_reff_profiles %>% 
  ggplot(aes(x = W, y = Reff, col = school)) +
    geom_line(size = 1.2) +
    theme_classic() +
    ylim(c(0,2)) +
    geom_hline(yintercept = 1, lty = 2) +
    scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-2,1,100),
                       limits = c(1e-4, 1)) +
    labs(x = "mean worm burden (W)",
         y = expression(italic(R[eff])))