> library(magrittr)
> library(dplyr)
> library(tidyr)
> library(purrr)
> library(readr)The trick here is to replace the concentration column by two columns minc and maxc where maxc = amount_anti_mg/kg/week. Whenever there is only 1 antibiotic, minc is equal to amount_anti_mg/kg/week as well but equal to 0 otherwise. Whenever there is no antibiotics at all, both minc and maxc are set to 0. For simplicity, we also replace missing value of anti_name by the character string no_anti:
> (feeds <- read_csv("Feed_Chapter4_data_10Dec.csv") %>%
+ mutate(minc = if_else(option > 0, 0, `amount_anti_mg/kg/week`),
+ maxc = `amount_anti_mg/kg/week`) %>%
+ mutate_at(vars(minc, maxc), replace_na, 0) %>%
+ mutate_at(vars(anti_name), replace_na, "no_anti") %>%
+ group_by(farmid, flockid, weekid, anti_name) %>%
+ summarise_at(vars(minc, maxc), sum) %>%
+ ungroup())
# A tibble: 6,767 x 6
farmid flockid weekid anti_name minc maxc
<chr> <chr> <chr> <chr> <dbl> <dbl>
1 75-001 75-001_1 75-001_1_1 enramycin 15.4 15.4
2 75-001 75-001_1 75-001_1_10 no_anti 0 0
3 75-001 75-001_1 75-001_1_11 no_anti 0 0
4 75-001 75-001_1 75-001_1_12 no_anti 0 0
5 75-001 75-001_1 75-001_1_13 no_anti 0 0
6 75-001 75-001_1 75-001_1_14 enramycin 3.46 3.46
7 75-001 75-001_1 75-001_1_15 no_anti 0 0
8 75-001 75-001_1 75-001_1_16 no_anti 0 0
9 75-001 75-001_1 75-001_1_17 no_anti 0 0
10 75-001 75-001_1 75-001_1_18 no_anti 0 0
# … with 6,757 more rowsNote: option == 1 means there is a label problem.
This being done, it’s now fairly easy to compute, the range of the mean concentration of, for example, colistin:
> feeds %>%
+ filter(anti_name %in% c("colistin", "no_anti")) %>%
+ summarise(minc = mean(minc), maxc = mean(maxc))
# A tibble: 1 x 2
minc maxc
<dbl> <dbl>
1 0.345 0.875Similarly, the range of the mean concentration of colistin per flock:
> feeds %>%
+ filter(anti_name %in% c("colistin", "no_anti")) %>%
+ group_by(flockid) %>%
+ summarise(minc = mean(minc), maxc = mean(maxc))
# A tibble: 259 x 3
flockid minc maxc
<chr> <dbl> <dbl>
1 75-001_1 0 0
2 75-001_2 0 0
3 75-001_3 0 0
4 75-001_4 0 0
5 75-001_5 0 0
6 75-001_6 0 0
7 75-001_7 0 0
8 75-002_1 0 0
9 75-002_2 0 0
10 75-002_3 0 0
# … with 249 more rowsAnd the range of the mean concentration of colistin per farm:
> feeds %>%
+ filter(anti_name %in% c("colistin", "no_anti")) %>%
+ group_by(farmid) %>%
+ summarise(minc = mean(minc), maxc = mean(maxc))
# A tibble: 113 x 3
farmid minc maxc
<chr> <dbl> <dbl>
1 75-001 0 0
2 75-002 0 0
3 75-004 0 0
4 75-005 0 0
5 75-006 0 0
6 75-008 0 0
7 75-009 0 0
8 75-010 0 0
9 75-011 0 0
10 75-012 0 0
# … with 103 more rowsThis is the list of all antibiotics:
> anti_list <- setdiff(unique(feeds$anti_name), "no_anti")The following function computes the ranges for all the antibiotics:
> compute_range <- function(x) {
+ anti_list %>%
+ map_dfr(~ filter(x, anti_name %in% c(.x, "no_anti")) %>%
+ summarise(minc = mean(minc), maxc = mean(maxc))) %>%
+ mutate(antibiotic = anti_list) %>%
+ select(antibiotic, everything())
+ }Let’s try it for all the weeks pooled together:
> compute_range(feeds)
# A tibble: 8 x 3
antibiotic minc maxc
<chr> <dbl> <dbl>
1 enramycin 2.19 2.42
2 chlortetracycline 6.65 6.72
3 bacitracin 3.44 5.54
4 virginiamycin 0.102 0.140
5 avilamycin 0 0.165
6 flavomycin 0 0.0763
7 colistin 0.345 0.875
8 Oxytetracyline 0 0.0879By flock:
> feeds %>%
+ group_by(flockid) %>%
+ group_modify(~ compute_range(.x)) %>%
+ ungroup()
# A tibble: 2,800 x 4
flockid antibiotic minc maxc
<chr> <chr> <dbl> <dbl>
1 75-001_1 enramycin 2.28 2.28
2 75-001_1 chlortetracycline 0 0
3 75-001_1 bacitracin 0 0
4 75-001_1 virginiamycin 0 0
5 75-001_1 avilamycin 0 0
6 75-001_1 flavomycin 0 0
7 75-001_1 colistin 0 0
8 75-001_1 Oxytetracyline 0 0
9 75-001_2 enramycin 2.86 2.86
10 75-001_2 chlortetracycline 0 0
# … with 2,790 more rowsBy farm:
> feeds %>%
+ group_by(farmid) %>%
+ group_modify(~ compute_range(.x)) %>%
+ ungroup()
# A tibble: 1,008 x 4
farmid antibiotic minc maxc
<chr> <chr> <dbl> <dbl>
1 75-001 enramycin 3.65 3.65
2 75-001 chlortetracycline 0 0
3 75-001 bacitracin 0 0
4 75-001 virginiamycin 0 0
5 75-001 avilamycin 0 0
6 75-001 flavomycin 0 0
7 75-001 colistin 0 0
8 75-001 Oxytetracyline 0 0
9 75-002 enramycin 4.20 4.20
10 75-002 chlortetracycline 0 0
# … with 998 more rowsIf now we are interested in frequencies of usage instead of concentration, it’s basically the same as above, replacing minc and maxc by presence / absence data. For all weeks pooled together:
> feeds %>%
+ mutate_at(vars(minc, maxc), as.logical) %>%
+ compute_range()
# A tibble: 8 x 3
antibiotic minc maxc
<chr> <dbl> <dbl>
1 enramycin 0.365 0.397
2 chlortetracycline 0.213 0.215
3 bacitracin 0.131 0.181
4 virginiamycin 0.0125 0.0309
5 avilamycin 0 0.0198
6 flavomycin 0 0.0198
7 colistin 0.00587 0.0163
8 Oxytetracyline 0 0.00265By flock:
> feeds %>%
+ mutate_at(vars(minc, maxc), as.logical) %>%
+ group_by(flockid) %>%
+ group_modify(~ compute_range(.x)) %>%
+ ungroup()
# A tibble: 2,800 x 4
flockid antibiotic minc maxc
<chr> <chr> <dbl> <dbl>
1 75-001_1 enramycin 0.296 0.296
2 75-001_1 chlortetracycline 0 0
3 75-001_1 bacitracin 0 0
4 75-001_1 virginiamycin 0 0
5 75-001_1 avilamycin 0 0
6 75-001_1 flavomycin 0 0
7 75-001_1 colistin 0 0
8 75-001_1 Oxytetracyline 0 0
9 75-001_2 enramycin 0.375 0.375
10 75-001_2 chlortetracycline 0 0
# … with 2,790 more rowsBy farm:
> feeds %>%
+ mutate_at(vars(minc, maxc), as.logical) %>%
+ group_by(farmid) %>%
+ group_modify(~ compute_range(.x)) %>%
+ ungroup()
# A tibble: 1,008 x 4
farmid antibiotic minc maxc
<chr> <chr> <dbl> <dbl>
1 75-001 enramycin 0.531 0.531
2 75-001 chlortetracycline 0 0
3 75-001 bacitracin 0 0
4 75-001 virginiamycin 0 0
5 75-001 avilamycin 0 0
6 75-001 flavomycin 0 0
7 75-001 colistin 0 0
8 75-001 Oxytetracyline 0 0
9 75-002 enramycin 0.699 0.699
10 75-002 chlortetracycline 0 0
# … with 998 more rows