analysis.Rmd

---
title: "AMR selective pressure"
csl: the-american-naturalist.csl
output:
  html_document:
    theme: cerulean
    toc: yes
  pdf_document:
    toc: yes
<!-- bibliography: references.bib -->
editor_options: 
  chunk_output_type: console
---

<!--
IMAGES:
Insert them with: ![alt text](image.png)
You can also resize them if needed: convert image.png -resize 50% image.png
If you want to center the image, go through HTML code:
<div style="text-align:center"><img src ="image.png"/></div>

REFERENCES:
For references: Put all the bibTeX references in the file "references.bib"
in the current folder and cite the references as @key or [@key] in the text.
Uncomment the bibliography field in the above header and put a "References"
title wherever you want to display the reference list.
-->

<style type="text/css">
.main-container {
  max-width: 1370px;
  margin-left: auto;
  margin-right: auto;
}
</style>

```{r general_options, include = FALSE}
knitr::knit_hooks$set(
  margin = function(before, options, envir) {
    if (before) par(mgp = c(1.5, .5, 0), bty = "n", plt = c(.11, .97, .135, 1))
    else NULL
  },
  prompt = function(before, options, envir) {
    options(prompt = if (options$engine %in% c("sh", "bash")) "$ " else "> ")
  })

knitr::opts_chunk$set(margin = TRUE, prompt = TRUE, comment = "",
                      collapse = TRUE, cache = FALSE, autodep = TRUE,
                      dev.args = list(pointsize = 11), fig.height = 3.5,
                      fig.width = 4.24725, fig.retina = 2, fig.align = "center")

options(width = 137)
```

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

```{r}
library(magrittr)
library(dplyr)
```

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

```{r}
#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)
}
```

```{r fig.width = 4.5, fig.height = 4, eval = FALSE}
vietnam %>% 
  mutate(amu = `Bodymass (1000s tonnes)` * `mg/Kg`) %>% 
  arrange(desc(amu)) %>% 
  skylineplot()
```

```{r}
new_data <- setNames(read.csv("new_data.txt"),
                     c("Species", "Bodymass (1000s tonnes)", "mg/Kg", "lower", "upper"))
```

```{r fig.width = 4.5, fig.height = 4}
new_data %>% 
  mutate(amu = `Bodymass (1000s tonnes)` * `mg/Kg`) %>% 
  arrange(desc(amu)) %>% 
  skylineplot()
```

```{r fig.width = 6, fig.height = 4}
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

```{r}
vn <- readxl::read_excel("AMU_VN_Europe.xlsx", "VN summary")
```

```{r}
eu <- readxl::read_excel("AMU_VN_Europe.xlsx", "Europe summary")
```


```{r}
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))
```

```{r}
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))
```


```{r}
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))
}
```

```{r fig.width = 8}
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)
```

```{r}
vn2 <- vn
vn2[1, 2] <- 262.8
vn2[2, 2] <- 245.6
vn2[2, 3] <- 11125.5
```

```{r fig.width = 8}
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)
```