civ-query-predict

library(eidith)
library(tidyverse)
library(rnaturalearth)
library(sf)
library(leaflet)
library(leaflet.extras)
# library(mapedit)
# library(mapview)
# library(ggmap)
library(lubridate)
library(patchwork)
library(cowplot)
# eidith::import_local_db("eha_with_malaysia")

# We wish to present bat only results , as supp mat tableswith the follwoing details :
#   -the number of caught animals, per species and per site
#   -number of samples collected and screened per species, per viral family (not only corona but all)
ed2_events()
conn <- eidith:::eidith_db()
DBI::dbListTables(conn)

# animals sampling
e <-  DBI::dbReadTable(conn, "events_2") %>%
  as_tibble() %>%
  filter(country == "Ivory Coast") %>%
  filter(concurrent_sampling_site != "Independent Site") %>%
  select(gains4_event_id, concurrent_sampling_site, site_name, site_latitude, site_longitude)

a <- DBI::dbReadTable(conn, "animals_2") %>%
  as_tibble() %>%
  filter(country == "Ivory Coast") %>%
  filter(taxa_group == "bats")

a <- inner_join(e, a)

n_distinct(a$gains4_sample_unit_id)

a %>%
  group_by(season) %>%
  count()

a %>%
  group_by(genus, species) %>%
  count() %>%
  arrange(-n)

a %>%
  filter(species != "NULL") %>%
  nrow()

a %>%
  group_by(site_name) %>%
  count()%>%
  arrange(-n)

a %>%
  group_by(age_class) %>%
  count()

a %>%
  group_by(genus) %>%
  count() %>%
  ungroup() %>%
  mutate(freq = formattable::percent(n / sum(n), digits = 1L)) %>%
  mutate(genus = fct_reorder(genus, n)) %>%
  mutate(label = paste0(n, " (", freq, ")")) %>%
  ggplot(aes(x = genus, y = n)) +
  geom_bar(stat = "identity", fill = "gray70") +
  geom_text(aes(label=label, y = n+1), hjust = "left", size = 5) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 160)) +
  labs(x = "", y = "") +
  coord_flip() +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 14, color = "black"),
        axis.text.x = element_blank(),
        axis.title.x = element_text(size = 14, hjust = 0, color = "black"),
        axis.ticks = element_blank(),
        panel.grid = element_blank())

# test results
t <- DBI::dbReadTable(conn, "tests_2") %>%
  as_tibble() %>%
  filter(country == "Ivory Coast") %>%
  inner_join(a %>% select(animal_id, genus, species))

pos <- t %>%
  filter(confirmation_result == "Positive") %>%
  left_join(a)

pos$season

write_csv(a, here::here("inst/queries/anne_20201125/civ_bat_collection_data.csv"))
write_csv(t, here::here("inst/queries/anne_20201125/civ_bat_testing_data.csv"))

# humans

h <- DBI::dbReadTable(conn, "human_2") %>%
  as_tibble() %>%
  filter(country == "Ivory Coast") %>%
  inner_join(e)

no_bats <- h %>%
  filter(bats_contact == "none") %>%
  nrow()

(nrow(h) - no_bats) / nrow(h)

hunt_bats <- h %>%
  filter(str_detect(bats_contact, "hunted")) %>%
  nrow()

hunt_bats/(nrow(h) - no_bats)

s <- DBI::dbReadTable(conn, "specimens_2") %>%
  as_tibble() %>%
  filter(country == "Ivory Coast") %>%
  select(-starts_with("season")) %>%
  inner_join(h,  by = c("project", "country", "gains4_sample_unit_id"))

t <- DBI::dbReadTable(conn, "tests_2") %>%
  as_tibble() %>%
  filter(country == "Ivory Coast") %>%
  inner_join(s, by = c("project", "country", "gains4_specimen_id"))

n_distinct(h$gains4_sample_unit_id)
n_distinct(t$gains4_sample_unit_id)

h %>%
  group_by(gender) %>%
  count()

h %>%
  group_by(site_name, gender) %>%
  summarize(average_age = mean(as.numeric(age))) %>%
  ungroup() %>%
  pivot_wider(names_from = gender, values_from = average_age) %>%
  write_csv(here::here("inst/queries/anne_20201125/civ_human_age_by_gender_and_site.csv"))


h %>%
  group_by(primary_livelihood) %>%
   count(sort = TRUE) %>%
  write_csv(here::here("inst/queries/anne_20201125/civ_human_primary_livelihood_counts.csv"))

t %>%
  distinct(gains4_sample_unit_id, gender) %>%
  group_by(gender) %>%
  count()

write_csv(h, here::here("inst/queries/anne_20201125/civ_human_surveillance_data.csv"))
write_csv(t, here::here("inst/queries/anne_20201125/civ_human_testing_data.csv"))

##### unused mapping code

# locs <- a %>%
#   distinct(site_name, site_latitude, site_longitude) %>%
#   mutate_at(.vars = c("site_latitude", "site_longitude"), as.numeric)

# civ <- ne_countries(country = "Ivory Coast", scale = 'large') %>%
#   st_as_sf()
# africa <- ne_countries(continent = "Africa", scale = 'large') %>%
#   st_as_sf()

# leaflet() %>%
#   setView(zoom = 15, lng = mean(locs$site_longitude), lat = mean(locs$site_latitude))  %>%
#   addPolygons(data = africa, stroke = TRUE, weight = 0.75 ,color = "black",
#               fill = TRUE, fillColor = '#798188', fillOpacity = 0.75) %>%
#   addPolygons(data = civ, stroke = TRUE, weight = 0.75 ,color = "black",
#               fill = TRUE, fillColor = "#FFFF99", fillOpacity = 0.85,
#               label = "", labelOptions = labelOptions(noHide = TRUE,
#                                                                 textOnly = TRUE,
#                                                                 offset = c(0, -90),
#                                                                 textsize = "18px",
#                                                                 opacity = 1)) %>%
#   addCircles(data = locs,
#              lng = ~site_longitude, lat = ~site_latitude, label = ~site_name,
#              color = "black", opacity = 1,
#              labelOptions = labelOptions(noHide = T, textsize = "12px", textOnly = TRUE, offset = c(5,0))) %>%
#   addScaleBar() %>%
#   setMapWidgetStyle(list(background= "#CCE2F5"))

# pd <- a %>%
#   group_by(genus) %>%
#   mutate(n = n()) %>%
#   ungroup() %>%
#   arrange(-n) %>%
#   mutate(genus = fct_reorder(genus, n)) %>%
#   group_by(genus, site_name) %>%
#   count() %>%
#   ungroup()
#
# plots <- map(unique(pd$site_name), function(x){
#   ggplot(pd %>% filter(site_name == x), aes(x = genus, y = n)) +
#     geom_bar(stat = 'identity') +
#     coord_flip() +
#     theme_minimal()
# })