AMR selective pressure

Skyline plot
VN and EU summaries

AMU / kg of animal

AMU / kg of human

We need

AMU in human
AMU in animals
number of people
weight of people
number of animal heads
weight of animal

Total human biomass:

$m_h = \sum_a n_{h_a} \times m_{h_a}$

Total biomass of animal $a$ :

$m_c = \sum_a n_{c_a} \times m_{c_a}$

$AMU = DDD \times d \times p$

were

$DDD$ = Defined Daily Dose (in g / kg / day)
$d$ = duration of treatment (in days)
$p$ = probability of treatment in an individual’s lifetime (or number of treaments divided by number of individuals)

Skyline plot

> library(magrittr)
> library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

> vietnam <- readxl::read_excel("AMU_VN_Europe_2.xlsx")
> #vietnam$Species %<>% sapply(sub, pattern = "Pigs", replacement = "Pig")

> #colors <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a','#ffff99')
> colors <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#FFF933','#808080','#cab2d6','#6a3d9a','#ffff99')
> skylineplot <- function(df) {
+   library(magrittr)
+   df %$% 
+     structure(list(breaks = cumsum(c(0, `Bodymass (1000s tonnes)`)),
+                    counts = `mg/Kg`,
+                    equidist = TRUE),
+               class = "histogram") %>%
+     assign("tmp", ., 1) %>% 
+     plot(col = colors, xlab = "cumulative biomass (1,000 tonnes)                          ", axes = FALSE,
+          ylab = "AMU (mg/kg)", main = NA, xlim = c(0, 20000), ylim = c(0, 500))
+   ats <- seq(0, 15000, 5000)
+   axis(1, ats, sub("000$", ",000", ats))
+   axis(2)
+   legend("topright", legend = df$Species, fill = colors, bty = "n")
+   
+   mids <- tmp$breaks[-length(tmp$breaks)] + diff(tmp$breaks) / 2
+   arrows(mids, df$lower, mids, df$upper, .05, 90, 3)
+ }

> vietnam %>% 
+   mutate(amu = `Bodymass (1000s tonnes)` * `mg/Kg`) %>% 
+   arrange(desc(amu)) %>% 
+   skylineplot()

> new_data <- setNames(read.csv("new_data.txt"),
+                      c("Species", "Bodymass (1000s tonnes)", "mg/Kg", "lower", "upper"))

> new_data %>% 
+   mutate(amu = `Bodymass (1000s tonnes)` * `mg/Kg`) %>% 
+   arrange(desc(amu)) %>% 
+   skylineplot()

> new_data2 <- new_data
> new_data2[, 1] <- as.character(new_data2[, 1])
> 
> new_data2[8, 2] <- new_data2[8, 2] + new_data2[10, 2]
> new_data2 <- new_data2[-10, ]
> new_data2[8, 1] <- "Goat and sheep"
> 
> new_data2[7, 2] <- new_data2[7, 2] + new_data2[9, 2]
> new_data2[7, 1] <- "Other avian"
> new_data2 <- new_data2[-9, ]
> 
> new_data2[2, 1] <- "Pig"
> 
> new_data2[4, 3] <- mean(unlist(new_data2[4, 4:5]))
> new_data2[6, 3] <- mean(unlist(new_data2[6, 4:5]))
> new_data2[8, 3] <- mean(unlist(new_data2[8, 4:5]))
> 
> new_data2 %>% 
+   mutate(amu = `Bodymass (1000s tonnes)` * `mg/Kg`) %>% 
+   arrange(desc(amu)) %>% 
+   skylineplot()

VN and EU summaries

> vn <- readxl::read_excel("AMU_VN_Europe.xlsx", "VN summary")

> eu <- readxl::read_excel("AMU_VN_Europe.xlsx", "Europe summary")

> vn %>% 
+   mutate(amu = `mg/Kg (unadjusted)` * `Bodymass (1000s tonnes)`,
+          breaks = 100 * `Bodymass (1000s tonnes)` / sum(`Bodymass (1000s tonnes)`)) %>% 
+   arrange(desc(amu)) %$%
+   structure(list(breaks = cumsum(c(0, breaks)),
+                  counts = `mg/Kg (unadjusted)`,
+                  equidist = TRUE),
+             class = "histogram") %>% 
+   plot(col = colors, xlab = "proportion of total biomass", axes = T,
+           ylab = "AMU (mg/kg)", main = NA, ylim = c(0, 210))

> eu %>% 
+   mutate(amu = `mg/Kg (unadjusted)` * `Bodymass (1000s tonnes)`,
+          breaks = 100 * `Bodymass (1000s tonnes)` / sum(`Bodymass (1000s tonnes)`)) %>% 
+   arrange(desc(amu)) %$%
+   structure(list(breaks = cumsum(c(0, breaks)),
+                  counts = `mg/Kg (unadjusted)`,
+                  equidist = TRUE),
+             class = "histogram") %>% 
+   plot(col = colors, xlab = "proportion of total biomass", axes = T,
+           ylab = "AMU (mg/kg)", main = NA, ylim = c(0, 210))

> comparison_plot <- function(x) {
+   x %>% 
+     mutate(amu = `mg/Kg (unadjusted)` * `Bodymass (1000s tonnes)`,
+            breaks = 100 * `Bodymass (1000s tonnes)` / sum(`Bodymass (1000s tonnes)`)) %>% 
+     arrange(desc(amu)) %$%
+     structure(list(breaks = cumsum(c(0, breaks)),
+                    counts = `mg/Kg (unadjusted)`,
+                    equidist = TRUE),
+               class = "histogram") %>% 
+     plot(col = colors, xlab = "proportion of total biomass", axes = T,
+              ylab = "AMU (mg/kg)", main = NA, ylim = c(0, 300))
+ }

> opar <- par(mfrow = c(1, 2))
> comparison_plot(vn)
> mtext("Vietnam", 3, -1)
> text(36.33971, 100, "animals")
> text(86.33971, 100, "humans", col = "white")
> comparison_plot(eu)
> mtext("European Union", 3, -1)
> text(31.83558, 100, "animals")
> text(81.83558, 100, "humans", col = "white")

> par(opar)

> vn2 <- vn
> vn2[1, 2] <- 262.8
> vn2[2, 2] <- 245.6
> vn2[2, 3] <- 11125.5

> opar <- par(mfrow = c(1, 2), plt = c(.11, .97, .135, .98))
> comparison_plot(vn2)
> mtext("Vietnam", 3, -1)
> text(36.33971, 100, "animals")
> text(86.33971, 100, "humans", col = "white")
> comparison_plot(eu)
> mtext("European Union", 3, -1)
> text(31.83558, 100, "animals")
> text(81.83558, 100, "humans", col = "white")

> par(opar)