code to make heatmap of CoV prevalence by bat family and quarter of t…

…he year

R/plot_quarterly_heatmap.R

@@ -0,0 +1,33 @@
+#' Make a heatmap of CoV prevalence by bat family and quarter of the year
+#'
+#'
+#' @title plot_quarterly_heatmap
+#' @param heat_data summary table with CoV prevalence by bat family and quarter of the year
+#' @return 
+#' @author Cecilia Sanchez
+
+plot_quarterly_heatmap <- function(heat_data){
+
+  ggplot(data = heat_data, mapping = aes(x = quarter, y = family,
+                                         fill = prop)) +
+    geom_tile(color = "black") +
+    scale_fill_viridis_c(name = "Prop. \nCoV +", option = "magma",
+                         direction = -1) +
+    # have text be labeled white if background fill is dark
+    geom_text(aes(label = paste0(n, "/", totN), color = prop > 0.15), size = 5) +
+    scale_color_manual(guide = "none", values = c("black", "white")) +
+    scale_x_discrete(labels = c("Jan-Mar", "Apr-Jun", "Jul-Sep", "Oct-Dec")) +
+    ylab("") + xlab("") +
+    theme_bw() +
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          axis.text = element_text(color = "black", size = 14),
+          legend.text = element_text(size = 14),
+          legend.title = element_text(size = 16),
+          legend.box.spacing = unit(0, "pt"),
+          legend.margin = margin(0,5,0,5),
+          plot.margin = margin(t = 2, r = 0, b = -10, l = -10))
+
+  ggsave("outputs/quarterly_prevalence_heatmap.png", width = 6, height = 4,
+         units = "in", dpi = 600)  
+}

R/prep_quarterly_heatmap_data.R

@@ -0,0 +1,43 @@
+#' Prep data to create a heatmap of CoV prevalence by bat family and quarter of the year
+#'
+#'
+#' @title prep_quarterly_heatmap_data
+#' @param CoV_prev dataset containing info on capture date, bat family, and CoV test result
+#' @return 
+#' @author Cecilia Sanchez
+
+prep_quarterly_heatmap_data <- function(CoV_prev){
+
+  # do some initial data cleaning
+  cleaned <- CoV_prev %>% 
+    janitor::clean_names() %>% 
+    mutate(date = lubridate::dmy(date),
+           quarter = as.factor(lubridate::quarter(date))) %>% 
+    rename(cov_result = co_v_final_results) %>% 
+    filter(!is.na(cov_result)) %>% 
+    mutate(cov_result = as.factor(cov_result))
+
+  # summarize proportion positive for each family and quarter
+  quarter_family_summary <- cleaned %>% 
+    group_by(family, quarter, cov_result, .drop = F) %>%
+    dplyr::summarize(n = n()) %>%
+    mutate(prop = n / sum(n)) %>%
+    mutate(totN = sum(n)) %>% 
+    filter(cov_result == "Positive")
+
+  # also calculate a summary only by quarter, lumping all families together
+  quarter_summary <- cleaned %>% 
+    group_by(quarter, cov_result, .drop = F) %>%
+    dplyr::summarize(n = n()) %>%
+    mutate(prop = n / sum(n)) %>%
+    mutate(totN = sum(n)) %>% 
+    filter(cov_result == "Positive") %>% 
+    mutate(family = "Overall")
+
+  # bind everything together to plot
+  heat_data <- bind_rows(quarter_family_summary, quarter_summary) %>%
+    mutate_at(vars(family), as.factor) %>%
+    mutate(family = fct_relevel(family, "Overall")) %>% 
+    drop_na()
+
+}

_targets.R

@@ -41,7 +41,10 @@ data_input_targets <- tar_plan(
 
   CoV_species = read.csv(CoV_species_file, na.strings = c("NA", "n/a")),
   WABNet_coords = read.csv(WABNet_coords_file, na.strings = c("NA", "n/a")),
-  darkcides = read.csv(darkcides_file, na.strings = c("N/A"))
+  darkcides = read.csv(darkcides_file, na.strings = c("N/A")),
+
+  tar_file(CoV_prev_file, "data/WABNet_CoV_prevalence_01December2023.csv"),
+  CoV_prev = read.csv(CoV_prev_file, na.strings = "")
 )
 
 
@@ -70,7 +73,9 @@ data_processing_targets <- tar_plan(
     filter(name %in% species_for_enm),
 
   iucn_ranges = terra::wrap(get_iucn(IUCN_data = mammals_file, 
-                                     species = species_for_enm))
+                                     species = species_for_enm)),
+
+  heat_data = prep_quarterly_heatmap_data(CoV_prev)
 
 )
 
@@ -107,7 +112,8 @@ analysis_targets <- tar_plan(
 ## Outputs
 outputs_targets <- tar_plan(
 
-  multipanel = map_predicted_distributions()
+  multipanel = map_predicted_distributions(),
+  quarterly_heatmap = plot_quarterly_heatmap(heat_data)
 )
 
 

packages.R

@@ -30,6 +30,7 @@ library(conflicted)
 # library(dotenv)
 # library(usethis)
 # library(kableExtra)
+library(janitor)
 library(qs)
 library(spocc)
 library(CoordinateCleaner)