AMR/vignettes/AMR.Rmd

---
title: "How to conduct AMR data analysis"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{How to conduct AMR data analysis}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r setup, include = FALSE, results = 'markup'}
knitr::opts_chunk$set(
  warning = FALSE,
  collapse = TRUE,
  comment = "#>",
  fig.width = 7.5,
  fig.height = 5
)
```

**Note:** values on this page will change with every website update since they are based on randomly created values and the page was written in [R Markdown](https://rmarkdown.rstudio.com/). However, the methodology remains unchanged. This page was generated on `r format(Sys.Date(), "%d %B %Y")`.

# Introduction

Conducting AMR data analysis unfortunately requires in-depth knowledge from different scientific fields, which makes it hard to do right. At least, it requires:

* Good questions (always start with those!) and reliable data
* A thorough understanding of (clinical) epidemiology, to understand the clinical and epidemiological relevance and possible bias of results
* A thorough understanding of (clinical) microbiology/infectious diseases, to understand which microorganisms are causal to which infections and the implications of pharmaceutical treatment, as well as understanding intrinsic and acquired microbial resistance
* Experience with data analysis with microbiological tests and their results, to understand the determination and limitations of MIC values and their interpretations to SIR values
* Availability of the biological taxonomy of microorganisms and probably normalisation factors for pharmaceuticals, such as defined daily doses (DDD)
* Available (inter-)national guidelines, and profound methods to apply them

Of course, we cannot instantly provide you with knowledge and experience. But with this `AMR` package, we aimed at providing (1) tools to simplify antimicrobial resistance data cleaning, transformation and analysis, (2) methods to easily incorporate international guidelines and (3) scientifically reliable reference data, including the requirements mentioned above.

The `AMR` package enables standardised and reproducible AMR data analysis, with the application of evidence-based rules, determination of first isolates, translation of various codes for microorganisms and antimicrobial agents, determination of (multi-drug) resistant microorganisms, and calculation of antimicrobial resistance, prevalence and future trends.

# Preparation 

For this tutorial, we will create fake demonstration data to work with. 

You can skip to [Cleaning the data](#cleaning-the-data) if you already have your own data ready. If you start your analysis, try to make the structure of your data generally look like this:

```{r example table, echo = FALSE, results = 'asis'}
knitr::kable(
  data.frame(
    date = Sys.Date(),
    patient_id = c("abcd", "abcd", "efgh"),
    mo = "Escherichia coli",
    AMX = c("S", "S", "R"),
    CIP = c("S", "R", "S"),
    stringsAsFactors = FALSE
  ),
  align = "c"
)
``` 

## Needed R packages

As with many uses in R, we need some additional packages for AMR data analysis. Our package works closely together with the [tidyverse packages](https://www.tidyverse.org) [`dplyr`](https://dplyr.tidyverse.org/) and [`ggplot2`](https://ggplot2.tidyverse.org) by RStudio. The tidyverse tremendously improves the way we conduct data science - it allows for a very natural way of writing syntaxes and creating beautiful plots in R.

We will also use the `cleaner` package, that can be used for cleaning data and creating frequency tables.

```{r lib packages, message = FALSE, warning = FALSE, results = 'asis'}
library(dplyr)
library(ggplot2)
library(AMR)

# (if not yet installed, install with:)
# install.packages(c("dplyr", "ggplot2", "AMR"))
```

The `AMR` package contains a data set `example_isolates_unclean`, which might look data that users have extracted from their laboratory systems:

```{r}
example_isolates_unclean

# we will use 'our_data' as the data set name for this tutorial
our_data <- example_isolates_unclean
```

For AMR data analysis, we would like the microorganism column to contain valid, up-to-date taxonomy, and the antibiotic columns to be cleaned as SIR values as well.

## Taxonomy of microorganisms

With `as.mo()`, users can transform arbitrary microorganism names or codes to current taxonomy. The `AMR` package contains up-to-date taxonomic data. To be specific, currently included data were retrieved on `r format(AMR:::TAXONOMY_VERSION$LPSN$accessed_date, "%d %b %Y")`. 

The codes of the AMR packages that come from `as.mo()` are short, but still human readable. More importantly, `as.mo()`  supports all kinds of input:

```{r, message = FALSE}
as.mo("Klebsiella pneumoniae")
as.mo("K. pneumoniae")
as.mo("KLEPNE")
as.mo("KLPN")
```

The first character in above codes denote their taxonomic kingdom, such as Bacteria (B), Fungi (F), and Protozoa (P).

The `AMR` package also contain functions to directly retrieve taxonomic properties, such as the name, genus, species, family, order, and even Gram-stain. They all start with `mo_` and they use `as.mo()` internally, so that still any arbitrary user input can be used:

```{r, message = FALSE}
mo_family("K. pneumoniae")
mo_genus("K. pneumoniae")
mo_species("K. pneumoniae")

mo_gramstain("Klebsiella pneumoniae")

mo_ref("K. pneumoniae")

mo_snomed("K. pneumoniae")
```

Now we can thus clean our data:

```{r, echo = FALSE, message = FALSE}
mo_reset_session()
```

```{r, message = TRUE}
our_data$bacteria <- as.mo(our_data$bacteria, info = TRUE)
```

Apparently, there was some uncertainty about the translation to taxonomic codes. Let's check this:

```{r}
mo_uncertainties()
```

That's all good.

## Antibiotic results

The column with antibiotic test results must also be cleaned. The `AMR` package comes with three new data types to work with such test results: `mic` for minimal inhibitory concentrations (MIC), `disk` for disk diffusion diameters, and `sir` for SIR data that have been interpreted already. This package can also determine SIR values based on MIC or disk diffusion values, read more about that on the `as.sir()` page.

For now, we will just clean the SIR columns in our data using dplyr:

```{r}
# method 1, be explicit about the columns:
our_data <- our_data %>%
  mutate_at(vars(AMX:GEN), as.sir)

# method 2, let the AMR package determine the eligible columns
our_data <- our_data %>%
  mutate_if(is_sir_eligible, as.sir)

# result:
our_data
```

This is basically it for the cleaning, time to start the data inclusion.

## First isolates

We need to know which isolates we can *actually* use for analysis without repetition bias.

To conduct an analysis of antimicrobial resistance, you must [only include the first isolate of every patient per episode](https:/pubmed.ncbi.nlm.nih.gov/17304462/) (Hindler *et al.*, Clin Infect Dis. 2007). If you would not do this, you could easily get an overestimate or underestimate of the resistance of an antibiotic. Imagine that a patient was admitted with an MRSA and that it was found in 5 different blood cultures the following weeks (yes, some countries like the Netherlands have these blood drawing policies). The resistance percentage of oxacillin of all \emph{S. aureus} isolates would be overestimated, because you included this MRSA more than once. It would clearly be [selection bias](https://en.wikipedia.org/wiki/Selection_bias).

The Clinical and Laboratory Standards Institute (CLSI) appoints this as follows:

> *(...) When preparing a cumulative antibiogram to guide clinical decisions about empirical antimicrobial therapy of initial infections, **only the first isolate of a given species per patient, per analysis period (eg, one year) should be included, irrespective of body site, antimicrobial susceptibility profile, or other phenotypical characteristics (eg, biotype)**. The first isolate is easily identified, and cumulative antimicrobial susceptibility test data prepared using the first isolate are generally comparable to cumulative antimicrobial susceptibility test data calculated by other methods, providing duplicate isolates are excluded.*
<br>[M39-A4 Analysis and Presentation of Cumulative Antimicrobial Susceptibility Test Data, 4th Edition. CLSI, 2014. Chapter 6.4](https://clsi.org/standards/products/microbiology/documents/m39/)

This `AMR` package includes this methodology with the `first_isolate()` function and is able to apply the four different methods as defined by [Hindler *et al.* in 2007](https://academic.oup.com/cid/article/44/6/867/364325): phenotype-based, episode-based, patient-based, isolate-based. The right method depends on your goals and analysis, but the default phenotype-based method is in any case the method to properly correct for most duplicate isolates. Read more about the methods on the `first_isolate()` page.

The outcome of the function can easily be added to our data:

```{r 1st isolate}
our_data <- our_data %>%
  mutate(first = first_isolate(info = TRUE))
```

So only `r round((sum(our_data$first) / nrow(our_data) * 100))`% is suitable for resistance analysis! We can now filter on it with the `filter()` function, also from the `dplyr` package:

```{r 1st isolate filter}
our_data_1st <- our_data %>%
  filter(first == TRUE)
```

For future use, the above two syntaxes can be shortened:

```{r 1st isolate filter 2}
our_data_1st <- our_data %>%
  filter_first_isolate()
```

So we end up with `r format(nrow(our_data_1st), big.mark = " ")` isolates for analysis. Now our data looks like:

```{r preview data set 3}
our_data_1st
```

Time for the analysis.

# Analysing the data

The base R `summary()` function gives a good first impression, as it comes with support for the new `mo` and `sir` classes that we now have in our data set:

```{r}
summary(our_data_1st)

glimpse(our_data_1st)

# number of unique values per column:
sapply(our_data_1st, n_distinct)
```

## Availability of species

To just get an idea how the species are distributed, create a frequency table with `count()` based on the name of the microorganisms:

```{r freq 1}
our_data %>%
  count(mo_name(bacteria), sort = TRUE)

our_data_1st %>%
  count(mo_name(bacteria), sort = TRUE)
```

## Select and filter with antibiotic selectors

Using so-called antibiotic class selectors, you can select or filter columns based on the antibiotic class that your antibiotic results are in:

```{r bug_drg 2a}
our_data_1st %>%
  select(date, aminoglycosides())

our_data_1st %>%
  select(bacteria, betalactams())

our_data_1st %>%
  select(bacteria, where(is.sir))

# filtering using AB selectors is also possible:
our_data_1st %>%
  filter(any(aminoglycosides() == "R"))

our_data_1st %>%
  filter(all(betalactams() == "R"))

# even works in base R (since R 3.0):
our_data_1st[all(betalactams() == "R"), ]
```

## Generate antibiograms

Since AMR v2.0 (March 2023), it is very easy to create different types of antibiograms, with support for 20 different languages.

There are four antibiogram types, as proposed by Klinker *et al.* (2021, [DOI 10.1177/20499361211011373](https://doi.org/10.1177/20499361211011373)), and they are all supported by the new `antibiogram()` function:

1. **Traditional Antibiogram (TA)** e.g, for the susceptibility of *Pseudomonas aeruginosa* to piperacillin/tazobactam (TZP)
2. **Combination Antibiogram (CA)** e.g, for the sdditional susceptibility of *Pseudomonas aeruginosa* to TZP + tobramycin versus TZP alone
3. **Syndromic Antibiogram (SA)** e.g, for the susceptibility of *Pseudomonas aeruginosa* to TZP among respiratory specimens (obtained among ICU patients only)
4. **Weighted-Incidence Syndromic Combination Antibiogram (WISCA)** e.g, for the susceptibility of *Pseudomonas aeruginosa* to TZP among respiratory specimens (obtained among ICU patients only) for male patients age >=65 years with heart failure

In this section, we show how to use the `antibiogram()` function to create any of the above antibiogram types. For starters, this is what the included `example_isolates` data set looks like:

```{r}
example_isolates
```

### Traditional Antibiogram

To create a traditional antibiogram, simply state which antibiotics should be used. The `antibiotics` argument in the `antibiogram()` function supports any (combination) of the previously mentioned antibiotic class selectors:

```{r trad}
antibiogram(example_isolates,
            antibiotics = c(aminoglycosides(), carbapenems()))
```

Notice that the `antibiogram()` function automatically prints in the right format when using Quarto or R Markdown (such as this page), and even applies italics for taxonomic names (by using `italicise_taxonomy()` internally).

It also uses the language of your OS if this is either `r AMR:::vector_or(vapply(FUN.VALUE = character(1), AMR:::LANGUAGES_SUPPORTED_NAMES, function(x) x$exonym), quotes = FALSE, sort = FALSE)`. In this next example, we force the language to be Spanish using the `language` argument:

```{r trad2}
antibiogram(example_isolates,
            mo_transform = "gramstain",
            antibiotics = aminoglycosides(),
            ab_transform = "name",
            language = "es")
```

### Combined Antibiogram

To create a combined antibiogram, use antibiotic codes or names with a plus `+` character like this:

```{r comb}
antibiogram(example_isolates,
            antibiotics = c("TZP", "TZP+TOB", "TZP+GEN"))
```

### Syndromic Antibiogram

To create a syndromic antibiogram, the `syndromic_group` argument must be used. This can be any column in the data, or e.g. an `ifelse()` with calculations based on certain columns:

```{r synd}
antibiogram(example_isolates,
            antibiotics = c(aminoglycosides(), carbapenems()),
            syndromic_group = "ward")
```

### Weighted-Incidence Syndromic Combination Antibiogram (WISCA)

To create a WISCA, you must state combination therapy in the `antibiotics` argument (similar to the Combination Antibiogram), define a syndromic group with the `syndromic_group` argument (similar to the Syndromic Antibiogram) in which cases are predefined based on clinical or demographic characteristics (e.g., endocarditis in 75+ females).  This next example is a simplification without clinical characteristics, but just gives an idea of how a WISCA can be created:

```{r wisca}
wisca <- antibiogram(example_isolates,
                     antibiotics = c("AMC", "AMC+CIP", "TZP", "TZP+TOB"),
                     mo_transform = "gramstain",
                     minimum = 10, # this should be >= 30, but now just as example
                     syndromic_group = ifelse(example_isolates$age >= 65 &
                                                example_isolates$gender == "M",
                                              "WISCA Group 1", "WISCA Group 2"))
wisca
```

### Plotting antibiograms

Antibiograms can be plotted using `autoplot()` from the `ggplot2` packages, since this `AMR` package provides an extension to that function:

```{r}
autoplot(wisca)
```

To calculate antimicrobial resistance in a more sensible way, also by correcting for too few results, we use the `resistance()` and `susceptibility()` functions.

## Resistance percentages

The functions `resistance()` and `susceptibility()` can be used to calculate antimicrobial resistance or susceptibility. For more specific analyses, the functions `proportion_S()`, `proportion_SI()`, `proportion_I()`, `proportion_IR()` and `proportion_R()` can be used to determine the proportion of a specific antimicrobial outcome. 

All these functions contain a `minimum` argument, denoting the minimum required number of test results for returning a value. These functions will otherwise return `NA`. The default is `minimum = 30`, following the [CLSI M39-A4 guideline](https://clsi.org/standards/products/microbiology/documents/m39/) for applying microbial epidemiology.

As per the EUCAST guideline of 2019, we calculate resistance as the proportion of R (`proportion_R()`, equal to `resistance()`) and susceptibility as the proportion of S and I (`proportion_SI()`, equal to `susceptibility()`). These functions can be used on their own:

```{r}
our_data_1st %>% resistance(AMX)
```

Or can be used in conjunction with `group_by()` and `summarise()`, both from the `dplyr` package:

```{r}
our_data_1st %>%
  group_by(hospital) %>%
  summarise(amoxicillin = resistance(AMX))
```

----

*Author: Dr. Matthijs Berends, 26th Feb 2023*