Chapter 5 EDM-Dengue

5.1 Init

library(tidyverse)
library(lubridate)
library(rEDM)
library(janitor)
library(purrr)
library(cowplot)
library(ggsci)
library(naniar)

5.2 Data

dat <- read.csv("./_data/LOR_dengue_clim.csv") %>%
  mutate(date = ymd(date))

vis_miss(dat %>%
           arrange(year, month))

5.3 EDA

dat_proc <- dat %>%
  select_if(~ !any(is.na(.)))

vis_miss(dat_proc %>%
           arrange(year, month))

top8 <- dat %>%
  group_by(district) %>%
  summarise(cases = sum(den_cases)) %>%
  arrange(-cases) %>%
  slice(1:8) %>%
  select(district) %>%
  unlist()

capitals <- c("IQUITOS","YURIMAGUAS", "NAUTA", "RAMON CASTILLA",
              "REQUENA", "CONTAMANA", "BARRANCA", "PUTUMAYO")

dat_proc %>%
  filter(district %in% top8) %>%
  gather(var, val, den_cases, pdsi:tmmx) %>%
  mutate(date = ymd(paste(year, month, "1", sep = "/"))) %>%
  ggplot(aes(x = date, y = val, col = var)) +
  geom_line() +
  guides(col = F) +
  facet_wrap(var~district,scales = "free", ncol = 8) +
  theme_bw()

5.4 EDM

vars <- c("pdsi", "pr", "soil", "tmmx")
e <- 1:8

dat_n <- dat_proc %>%
  filter(district %in% top8) %>%
  select(-date) %>%
  group_by(UBIGEO, department, province, district, hogar,vivienda) %>%
  nest() %>%
  crossing(vars) %>%
  crossing(e) %>%
  #filter(district == "IQUITOS", vars == "pdsi") %>% #testing
  mutate(ccm = map2(.x = data, .y = vars,
                         .f = ~CCM(dataFrame = .x, E = 3, 
                                   columns = "den_cases",
                                   target = .y, libSizes = "10 70 5", sample = 100)),
         ccm2 = pmap(.l = list(a = data, b = vars, c = e),
                     .f = ~CCM(dataFrame = ..1, E = ..3, 
                               columns = "den_cases",
                               target = ..2, libSizes = "10 70 5", sample = 100)),
         title = map2_chr(.x = district, .y = vars, .f = ~paste0(.x, ":", .y)),
         title2 = pmap_chr(.l = list(a = district, b = vars, c = e),
                           .f = ~paste0(..1, ":", ..2,":E=",..3)),
         ccm_plot = map2(.x = ccm2,
                         .y = title2,
                        .f = ~plot_ccm(dat = .x, t = .y)))

#saveRDS(dat_n, "./_data/dat_n.rds")

dat_n <- readRDS("./_data/dat_n.rds")

plot_grid(plotlist = dat_n$ccm_plot, ncol = 28) 
ggsave("./_out/edm_plots.png", width = 49, height = 30)

e %>%
  map(~ggsave_ccm(dat = dat_n, out = .x))

top8 %>%
  map(~ggsave_ccm(dat = dat_n, out = .x, d = T))