Sensitivity analysis for county data

Last updated: 2019-07-24

Purpose

This code explores the sensitivity of our insecticide analyses to uncertainties in the original data sources. The main issues explored are:

• Certainty associated with LD50 estimates
• Interpolated values for acreage and insecticide use

Libraries & functions

library(tidyverse)

── Attaching packages ───────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

✔ ggplot2 3.2.0     ✔ purrr   0.3.2
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   0.8.3     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0

── Conflicts ──────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

Load data

pest <- read.csv("../output_big/bee_tox_index_cty.csv") %>%
  filter(YEAR < 2013 & YEAR > 1996)

pest_source <- read.csv("../output_big/bee_tox_index_cty_source.csv") %>%
  filter(YEAR < 2013 & YEAR > 1996)

Graph insecticide trends by source of LD50 data

pest_source_long <- pest_source %>%
  filter(cat == "I") %>%
  mutate(uncertainty = ifelse(source %in% c("epa","epa/eu","eu"),"low",
                          ifelse(source %in% c("cat_med","rac_med"),"high",
                                 "medium"))) %>%
  mutate(uncertainty = fct_relevel(uncertainty, "high", "medium", "low")) %>%
  select(YEAR, fips, source, uncertainty, kg_low, ct_tox_bil_low, or_tox_bil_low) %>%
  gather(key = "var",
         value = "amount",
         kg_low:or_tox_bil_low) %>%
  group_by(YEAR, uncertainty, var) %>%
  summarise(amount = sum(amount))

pest_source_wide <- pest_source_long %>%
  spread(key = uncertainty,
         value = amount) %>%
  mutate(perc_low = (low/(high + medium + low))*100,
          perc_high = (high/(high + medium + low))*100) %>%
  arrange(var)
  
pest_source_long$var <- fct_relevel(pest_source_long$var, "kg_low", "ct_tox_bil_low", "or_tox_bil_low")

var_names <- list(
  'kg_low'="Weight (mil kg)",
  'ct_tox_bil_low'=expression('Contact toxic load (10'^15*' LD'[50]*'s)'),
  'or_tox_bil_low'=expression('Oral toxic load (10'^15*' LD'[50]*'s)')
)

var_labeller <- function(variable,value){
  return(var_names[value])
}

text_lab <- data.frame(
  label = c("a)", "b)", "c)"),
  var   = c("kg_low", "ct_tox_bil_low", "or_tox_bil_low")
)

ggplot(pest_source_long, aes(fill = uncertainty, x = YEAR, y = amount/10^6)) +
  geom_col() +
  facet_wrap(~var, scales = "free_y",
         strip.position = "left",
         labeller = var_labeller) +
  ylab(NULL) +
  xlab(NULL) +
  theme_classic(base_size=18) +
  theme(strip.background = element_blank(),
        strip.placement = "outside") +
  scale_fill_brewer(palette = "Set1") +
  guides(fill=guide_legend(title=NULL)) 

Warning: The labeller API has been updated. Labellers taking `variable`and
`value` arguments are now deprecated. See labellers documentation.

Pesticides interpolated

# create summary table
pest_interp_reg <- pest %>%
  group_by(YEAR, region_name, interp_pest) %>%
  summarise(ct_tox = sum(ct_tox_bil_low, na.rm=TRUE),
            or_tox = sum(or_tox_bil_low, na.rm=TRUE),
            kg_low = sum(kg_low, na.rm=TRUE),
            n_interp = n()) %>%
  na.omit() %>%
  gather(value = "value",
         key = "var",
         ct_tox:n_interp) %>%
  spread(value = value,
         key = interp_pest) %>%
  replace_na(replace = list(yes = 0)) %>%
  mutate(total = no + yes,
         perc_interp = (yes/total)*100) %>%
  arrange(var)

pest_interp_nat <- pest %>%
  group_by(YEAR, interp_pest) %>%
  summarise(ct_tox = sum(ct_tox_bil_low, na.rm=TRUE),
            or_tox = sum(or_tox_bil_low, na.rm=TRUE),
            kg_low = sum(kg_low, na.rm=TRUE),
            n_interp = n()) %>%
  na.omit() %>%
  gather(value = "value",
         key = "var",
         ct_tox:n_interp) %>%
  spread(value = value,
         key = interp_pest) %>%
  replace_na(replace = list(yes = 0)) %>%
  mutate(total = no + yes,
         perc_interp = (yes/total)*100,
         region_name = "All regions (U.S.)") %>%
  arrange(var)  

pest_interp_tot <- rbind(pest_interp_reg, pest_interp_nat)

Warning in bind_rows_(x, .id): binding factor and character vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

ggplot(data = pest_interp_tot, aes(x = YEAR, y = perc_interp, colour = var)) +
  geom_line() +
  facet_wrap(~region_name) +
  theme_classic() +
  ylab("Interpolated values (% of total)") +
  xlab("") +
  theme(legend.title = element_blank()) +
  ggtitle("Contribution of counties with interpolated insecticide values")

Area interpolated

# create summary table
ha_interp_reg <- pest %>%
  group_by(YEAR, region_name, crop_ha_interp) %>%
  summarise(crop_ha = sum(crop_ha, na.rm=TRUE),
            n_interp = n()) %>%
  na.omit() %>%
  gather(value = "value",
         key = "var",
         crop_ha:n_interp) %>%
  spread(value = value,
         key = crop_ha_interp) %>%
  replace_na(replace = list(yes = 0)) %>%
  mutate(total = no + yes,
         perc_interp = (yes/total)*100) %>%
  arrange(var, perc_interp)

Warning: Factor `crop_ha_interp` contains implicit NA, consider using
`forcats::fct_explicit_na`

ha_interp_nat <- pest %>%
  group_by(YEAR, crop_ha_interp) %>%
  summarise(crop_ha = sum(crop_ha, na.rm=TRUE),
            n_interp = n()) %>%
  na.omit() %>%
  gather(value = "value",
         key = "var",
         crop_ha:n_interp) %>%
  spread(value = value,
         key = crop_ha_interp) %>%
  replace_na(replace = list(yes = 0)) %>%
  mutate(total = no + yes,
         perc_interp = (yes/total)*100,
         region_name = "All states (U.S.)") %>%
  arrange(var, perc_interp)

Warning: Factor `crop_ha_interp` contains implicit NA, consider using
`forcats::fct_explicit_na`

ha_interp_tot <- rbind(ha_interp_reg, ha_interp_nat)

Warning in bind_rows_(x, .id): binding factor and character vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

ggplot(data = ha_interp_tot, aes(x = YEAR, y = perc_interp, colour = var)) +
  geom_line() +
  facet_wrap(~region_name) +
  theme_classic() +
  ylab("Interpolated values (% of total)") +
  xlab("") +
  theme(legend.title = element_blank()) +
  scale_x_continuous(limits = c(1996,2013), breaks = c(1997,2002,2007,2012)) +
  ggtitle("Contribution of counties with interpolated acreage values")

Area treated interpolated

# create summary table
trt_ha_interp_reg <- pest %>%
  group_by(YEAR, region_name, trt_ins_ha_interp) %>%
  summarise(trt_ins_ha = sum(trt_ins_ha, na.rm=TRUE),
            n_interp_trt = n()) %>%
  na.omit() %>%
  gather(value = "value",
         key = "var",
         trt_ins_ha:n_interp_trt) %>%
  spread(value = value,
         key = trt_ins_ha_interp) %>%
  replace_na(replace = list(yes = 0)) %>%
  mutate(total = no + yes,
         perc_interp = (yes/total)*100) %>%
  arrange(var, perc_interp)

Warning: Factor `trt_ins_ha_interp` contains implicit NA, consider using
`forcats::fct_explicit_na`

trt_ha_interp_nat <- pest %>%
  group_by(YEAR, trt_ins_ha_interp) %>%
  summarise(trt_ins_ha = sum(trt_ins_ha, na.rm=TRUE),
            n_interp_trt = n()) %>%
  na.omit() %>%
  gather(value = "value",
         key = "var",
         trt_ins_ha:n_interp_trt) %>%
  spread(value = value,
         key = trt_ins_ha_interp) %>%
  replace_na(replace = list(yes = 0)) %>%
  mutate(total = no + yes,
         perc_interp = (yes/total)*100,
         region_name = "All states (U.S.)") %>%
  arrange(var, perc_interp)

Warning: Factor `trt_ins_ha_interp` contains implicit NA, consider using
`forcats::fct_explicit_na`

trt_ha_interp_tot <- rbind(trt_ha_interp_reg, trt_ha_interp_nat)

Warning in bind_rows_(x, .id): binding factor and character vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

ggplot(data = trt_ha_interp_tot, aes(x = YEAR, y = perc_interp, colour = var)) +
  geom_line() +
  facet_wrap(~region_name) +
  theme_classic() +
  ylab("Interpolated values (% of total)") +
  xlab("") +
  theme(legend.title = element_blank()) +
  scale_x_continuous(limits = c(1996,2013), breaks = c(1997,2002,2007,2012)) +
  ggtitle("Contribution of counties with interpolated treated acreage values")

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.3     purrr_0.3.2    
[5] readr_1.3.1     tidyr_0.8.3     tibble_2.1.3    ggplot2_3.2.0  
[9] tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1         RColorBrewer_1.1-2 cellranger_1.1.0  
 [4] pillar_1.4.2       compiler_3.6.1     tools_3.6.1       
 [7] zeallot_0.1.0      digest_0.6.20      lubridate_1.7.4   
[10] jsonlite_1.6       evaluate_0.14      nlme_3.1-140      
[13] gtable_0.3.0       lattice_0.20-38    pkgconfig_2.0.2   
[16] rlang_0.4.0        cli_1.1.0          rstudioapi_0.10   
[19] yaml_2.2.0         haven_2.1.1        xfun_0.8          
[22] withr_2.1.2        xml2_1.2.0         httr_1.4.0        
[25] knitr_1.23         vctrs_0.2.0        hms_0.5.0         
[28] generics_0.0.2     fs_1.3.1           grid_3.6.1        
[31] tidyselect_0.2.5   glue_1.3.1         R6_2.4.0          
[34] readxl_1.3.1       rmarkdown_1.14     modelr_0.1.4      
[37] magrittr_1.5       backports_1.1.4    scales_1.0.0      
[40] htmltools_0.3.6    rvest_0.3.4        assertthat_0.2.1  
[43] colorspace_1.4-1   labeling_0.3       stringi_1.4.3     
[46] lazyeval_0.2.2     munsell_0.5.0      broom_0.5.2       
[49] crayon_1.3.4      

---

This R Markdown site was created with workflowr

---