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)

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