AMR/vignettes/resistance_predict.Rmd

132 lines
5.2 KiB
Plaintext
Raw Normal View History

2019-01-02 23:24:07 +01:00
---
title: "How to predict antimicrobial resistance"
author: "Matthijs S. Berends"
date: '`r format(Sys.Date(), "%d %B %Y")`'
2019-01-02 23:24:07 +01:00
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{How to predict antimicrobial resistance}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
editor_options:
chunk_output_type: console
---
```{r setup, include = FALSE, results = 'markup'}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#",
2019-02-09 22:16:24 +01:00
fig.width = 7.5,
2019-02-11 10:27:10 +01:00
fig.height = 4.75
2019-01-02 23:24:07 +01:00
)
```
2019-02-09 22:16:24 +01:00
## Needed R packages
2019-06-02 19:23:19 +02:00
As with many uses in R, we need some additional packages for AMR 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 Dr Hadley Wickham. 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.
2019-02-09 22:16:24 +01:00
Our `AMR` package depends on these packages and even extends their use and functions.
```{r lib packages, message = FALSE}
library(dplyr)
library(ggplot2)
library(AMR)
# (if not yet installed, install with:)
# install.packages(c("tidyverse", "AMR"))
```
## Prediction analysis
2019-02-11 10:27:10 +01:00
Our package contains a function `resistance_predict()`, which takes the same input as functions for [other AMR analysis](./AMR.html). Based on a date column, it calculates cases per year and uses a regression model to predict antimicrobial resistance.
2019-02-09 22:16:24 +01:00
It is basically as easy as:
```{r, eval = FALSE}
2019-05-10 16:44:59 +02:00
# resistance prediction of piperacillin/tazobactam (TZP):
resistance_predict(tbl = example_isolates, col_date = "date", col_ab = "TZP", model = "binomial")
2019-02-09 22:16:24 +01:00
# or:
example_isolates %>%
2019-08-08 15:52:07 +02:00
resistance_predict(col_ab = "TZP",
model "binomial")
2019-02-09 22:16:24 +01:00
2019-05-10 16:44:59 +02:00
# to bind it to object 'predict_TZP' for example:
predict_TZP <- example_isolates %>%
2019-08-08 15:52:07 +02:00
resistance_predict(col_ab = "TZP",
model = "binomial")
2019-02-09 22:16:24 +01:00
```
2019-02-11 10:27:10 +01:00
The function will look for a date column itself if `col_date` is not set.
When running any of these commands, a summary of the regression model will be printed unless using `resistance_predict(..., info = FALSE)`.
2019-02-09 22:16:24 +01:00
```{r, echo = FALSE}
predict_TZP <- example_isolates %>%
2019-08-08 15:52:07 +02:00
resistance_predict(col_ab = "TZP", model = "binomial")
2019-02-09 22:16:24 +01:00
```
2019-02-11 10:27:10 +01:00
This text is only a printed summary - the actual result (output) of the function is a `data.frame` containing for each year: the number of observations, the actual observed resistance, the estimated resistance and the standard error below and above the estimation:
2019-02-09 22:16:24 +01:00
```{r}
2019-05-10 16:44:59 +02:00
predict_TZP
2019-02-09 22:16:24 +01:00
```
The function `plot` is available in base R, and can be extended by other packages to depend the output based on the type of input. We extended its function to cope with resistance predictions:
2019-02-11 10:27:10 +01:00
```{r, fig.height = 5.5}
2019-05-10 16:44:59 +02:00
plot(predict_TZP)
2019-02-09 22:16:24 +01:00
```
2019-02-11 10:27:10 +01:00
This is the fastest way to plot the result. It automatically adds the right axes, error bars, titles, number of available observations and type of model.
We also support the `ggplot2` package with our custom function `ggplot_rsi_predict()` to create more appealing plots:
2019-02-09 22:16:24 +01:00
```{r}
2019-05-10 16:44:59 +02:00
ggplot_rsi_predict(predict_TZP)
2019-02-11 10:27:10 +01:00
# choose for error bars instead of a ribbon
2019-05-10 16:44:59 +02:00
ggplot_rsi_predict(predict_TZP, ribbon = FALSE)
2019-02-09 22:16:24 +01:00
```
### Choosing the right model
Resistance is not easily predicted; if we look at vancomycin resistance in Gram positives, the spread (i.e. standard error) is enormous:
```{r}
example_isolates %>%
2019-06-22 21:33:13 +02:00
filter(mo_gramstain(mo, language = NULL) == "Gram-positive") %>%
2019-08-08 15:52:07 +02:00
resistance_predict(col_ab = "VAN", year_min = 2010, info = FALSE, model = "binomial") %>%
2019-02-11 10:27:10 +01:00
ggplot_rsi_predict()
2019-02-09 22:16:24 +01:00
```
Vancomycin resistance could be 100% in ten years, but might also stay around 0%.
2019-08-08 15:52:07 +02:00
You can define the model with the `model` parameter. The model chosen above is a generalised linear regression model using a binomial distribution, assuming that a period of zero resistance was followed by a period of increasing resistance leading slowly to more and more resistance.
2019-02-09 22:16:24 +01:00
Valid values are:
| Input values | Function used by R | Type of model |
|----------------------------------------|-------------------------------|-----------------------------------------------------|
| `"binomial"` or `"binom"` or `"logit"` | `glm(..., family = binomial)` | Generalised linear model with binomial distribution |
| `"loglin"` or `"poisson"` | `glm(..., family = poisson)` | Generalised linear model with poisson distribution |
| `"lin"` or `"linear"` | `lm()` | Linear model |
2019-02-11 10:27:10 +01:00
For the vancomycin resistance in Gram positive bacteria, a linear model might be more appropriate since no (left half of a) binomial distribution is to be expected based on the observed years:
2019-02-09 22:16:24 +01:00
```{r}
example_isolates %>%
2019-06-22 21:33:13 +02:00
filter(mo_gramstain(mo, language = NULL) == "Gram-positive") %>%
2019-05-10 16:44:59 +02:00
resistance_predict(col_ab = "VAN", year_min = 2010, info = FALSE, model = "linear") %>%
2019-02-11 10:27:10 +01:00
ggplot_rsi_predict()
2019-02-09 22:16:24 +01:00
```
This seems more likely, doesn't it?
2019-02-11 10:27:10 +01:00
The model itself is also available from the object, as an `attribute`:
```{r}
2019-05-10 16:44:59 +02:00
model <- attributes(predict_TZP)$model
2019-02-11 10:27:10 +01:00
summary(model)$family
summary(model)$coefficients
```