AMR/R/resistance_predict.R

389 lines
17 KiB
R
Raw Normal View History

2018-08-10 15:01:05 +02:00
# ==================================================================== #
# TITLE #
2020-10-08 11:16:03 +02:00
# Antimicrobial Resistance (AMR) Analysis for R #
2018-08-10 15:01:05 +02:00
# #
2019-01-02 23:24:07 +01:00
# SOURCE #
2020-07-08 14:48:06 +02:00
# https://github.com/msberends/AMR #
2018-08-10 15:01:05 +02:00
# #
# LICENCE #
# (c) 2018-2020 Berends MS, Luz CF et al. #
2020-10-08 11:16:03 +02:00
# Developed at the University of Groningen, the Netherlands, in #
# collaboration with non-profit organisations Certe Medical #
# Diagnostics & Advice, and University Medical Center Groningen. #
2018-08-10 15:01:05 +02:00
# #
2019-01-02 23:24:07 +01:00
# This R package is free software; you can freely use and distribute #
# it for both personal and commercial purposes under the terms of the #
# GNU General Public License version 2.0 (GNU GPL-2), as published by #
# the Free Software Foundation. #
# We created this package for both routine data analysis and academic #
# research and it was publicly released in the hope that it will be #
# useful, but it comes WITHOUT ANY WARRANTY OR LIABILITY. #
2020-10-08 11:16:03 +02:00
# #
# Visit our website for the full manual and a complete tutorial about #
# how to conduct AMR analysis: https://msberends.github.io/AMR/ #
2018-08-10 15:01:05 +02:00
# ==================================================================== #
#' Predict antimicrobial resistance
#'
#' Create a prediction model to predict antimicrobial resistance for the next years on statistical solid ground. Standard errors (SE) will be returned as columns `se_min` and `se_max`. See *Examples* for a real live example.
#' @inheritSection lifecycle Maturing lifecycle
#' @param col_ab column name of `x` containing antimicrobial interpretations (`"R"`, `"I"` and `"S"`)
2019-01-15 12:45:24 +01:00
#' @param col_date column name of the date, will be used to calculate years if this column doesn't consist of years already, defaults to the first column of with a date class
#' @param year_min lowest year to use in the prediction model, dafaults to the lowest year in `col_date`
2019-01-12 19:31:30 +01:00
#' @param year_max highest year to use in the prediction model, defaults to 10 years after today
#' @param year_every unit of sequence between lowest year found in the data and `year_max`
2018-08-10 15:01:05 +02:00
#' @param minimum minimal amount of available isolates per year to include. Years containing less observations will be estimated by the model.
#' @param model the statistical model of choice. This could be a generalised linear regression model with binomial distribution (i.e. using `glm(..., family = binomial)``, assuming that a period of zero resistance was followed by a period of increasing resistance leading slowly to more and more resistance. See Details for all valid options.
#' @param I_as_S a logical to indicate whether values `I` should be treated as `S` (will otherwise be treated as `R`). The default, `TRUE`, follows the redefinition by EUCAST about the interpretion of I (increased exposure) in 2019, see section *Interpretation of S, I and R* below.
#' @param preserve_measurements a logical to indicate whether predictions of years that are actually available in the data should be overwritten by the original data. The standard errors of those years will be `NA`.
#' @param info a logical to indicate whether textual analysis should be printed with the name and [summary()] of the statistical model.
2019-01-15 16:38:54 +01:00
#' @param main title of the plot
2019-02-11 10:27:10 +01:00
#' @param ribbon a logical to indicate whether a ribbon should be shown (default) or error bars
2019-05-13 16:35:48 +02:00
#' @param ... parameters passed on to functions
2019-11-29 19:43:23 +01:00
#' @inheritSection as.rsi Interpretation of R and S/I
2019-05-13 16:35:48 +02:00
#' @inheritParams first_isolate
#' @inheritParams graphics::plot
#' @details Valid options for the statistical model (parameter `model`) are:
#' - `"binomial"` or `"binom"` or `"logit"`: a generalised linear regression model with binomial distribution
#' - `"loglin"` or `"poisson"`: a generalised log-linear regression model with poisson distribution
#' - `"lin"` or `"linear"`: a linear regression model
#' @return A [data.frame] with extra class [`resistance_predict`] with columns:
#' - `year`
#' - `value`, the same as `estimated` when `preserve_measurements = FALSE`, and a combination of `observed` and `estimated` otherwise
#' - `se_min`, the lower bound of the standard error with a minimum of `0` (so the standard error will never go below 0%)
#' - `se_max` the upper bound of the standard error with a maximum of `1` (so the standard error will never go above 100%)
#' - `observations`, the total number of available observations in that year, i.e. \eqn{S + I + R}
#' - `observed`, the original observed resistant percentages
#' - `estimated`, the estimated resistant percentages, calculated by the model
#'
#' Furthermore, the model itself is available as an attribute: `attributes(x)$model`, please see *Examples*.
#' @seealso The [proportion()] functions to calculate resistance
#'
#' Models: [lm()] [glm()]
2018-08-10 15:01:05 +02:00
#' @rdname resistance_predict
#' @export
#' @importFrom stats predict glm lm
2019-01-02 23:24:07 +01:00
#' @inheritSection AMR Read more on our website!
2018-08-10 15:01:05 +02:00
#' @examples
#' x <- resistance_predict(example_isolates,
#' col_ab = "AMX",
#' year_min = 2010,
#' model = "binomial")
2019-01-15 12:45:24 +01:00
#' plot(x)
2020-05-16 21:40:50 +02:00
#' if (require("ggplot2")) {
#' ggplot_rsi_predict(x)
#' }
2018-08-10 15:01:05 +02:00
#'
2020-05-16 13:05:47 +02:00
#' # using dplyr:
2020-05-16 21:40:50 +02:00
#' if (require("dplyr")) {
2020-05-16 13:05:47 +02:00
#' x <- example_isolates %>%
#' filter_first_isolate() %>%
#' filter(mo_genus(mo) == "Staphylococcus") %>%
#' resistance_predict("PEN", model = "binomial")
#' plot(x)
2019-02-11 10:27:10 +01:00
#'
2020-05-16 13:05:47 +02:00
#' # get the model from the object
#' mymodel <- attributes(x)$model
#' summary(mymodel)
#' }
2019-02-11 10:27:10 +01:00
#'
#' # create nice plots with ggplot2 yourself
#' if (require("dplyr") & require("ggplot2")) {
#'
#' data <- example_isolates %>%
2018-12-22 22:39:34 +01:00
#' filter(mo == as.mo("E. coli")) %>%
2019-05-10 16:44:59 +02:00
#' resistance_predict(col_ab = "AMX",
2018-12-22 22:39:34 +01:00
#' col_date = "date",
2019-08-07 15:58:32 +02:00
#' model = "binomial",
2018-12-22 22:39:34 +01:00
#' info = FALSE,
2019-01-15 12:45:24 +01:00
#' minimum = 15)
2018-08-10 15:01:05 +02:00
#'
#' ggplot(data,
#' aes(x = year)) +
#' geom_col(aes(y = value),
#' fill = "grey75") +
#' geom_errorbar(aes(ymin = se_min,
#' ymax = se_max),
#' colour = "grey50") +
#' scale_y_continuous(limits = c(0, 1),
#' breaks = seq(0, 1, 0.1),
#' labels = paste0(seq(0, 100, 10), "%")) +
#' labs(title = expression(paste("Forecast of Amoxicillin Resistance in ",
2018-08-10 15:01:05 +02:00
#' italic("E. coli"))),
#' y = "%R",
2018-08-10 15:01:05 +02:00
#' x = "Year") +
#' theme_minimal(base_size = 13)
#' }
2019-05-13 16:35:48 +02:00
resistance_predict <- function(x,
2018-08-10 15:01:05 +02:00
col_ab,
2019-01-15 12:45:24 +01:00
col_date = NULL,
2018-08-10 15:01:05 +02:00
year_min = NULL,
year_max = NULL,
year_every = 1,
minimum = 30,
2019-08-07 15:37:39 +02:00
model = NULL,
2019-05-13 16:35:48 +02:00
I_as_S = TRUE,
2018-08-10 15:01:05 +02:00
preserve_measurements = TRUE,
info = interactive(),
2019-05-13 16:35:48 +02:00
...) {
meet_criteria(x, allow_class = "data.frame")
meet_criteria(col_ab, allow_class = "character", has_length = 1, is_in = colnames(x))
meet_criteria(col_date, allow_class = "character", has_length = 1, is_in = colnames(x), allow_NULL = TRUE)
meet_criteria(year_min, allow_class = c("numeric", "integer"), has_length = 1, allow_NULL = TRUE)
meet_criteria(year_max, allow_class = c("numeric", "integer"), has_length = 1, allow_NULL = TRUE)
meet_criteria(year_every, allow_class = c("numeric", "integer"), has_length = 1)
meet_criteria(minimum, allow_class = c("numeric", "integer"), has_length = 1)
meet_criteria(model, allow_class = c("character", "function"), has_length = 1, allow_NULL = TRUE)
meet_criteria(I_as_S, allow_class = "logical", has_length = 1)
meet_criteria(preserve_measurements, allow_class = "logical", has_length = 1)
meet_criteria(info, allow_class = "logical", has_length = 1)
2020-07-13 09:17:24 +02:00
stop_if(is.null(model), 'choose a regression model with the `model` parameter, e.g. resistance_predict(..., model = "binomial")')
2019-05-13 16:35:48 +02:00
dots <- unlist(list(...))
if (length(dots) != 0) {
# backwards compatibility with old parameters
dots.names <- dots %pm>% names()
2019-10-11 17:21:02 +02:00
if ("tbl" %in% dots.names) {
x <- dots[which(dots.names == "tbl")]
2019-05-13 16:35:48 +02:00
}
2019-10-11 17:21:02 +02:00
if ("I_as_R" %in% dots.names) {
2019-05-13 16:35:48 +02:00
warning("`I_as_R is deprecated - use I_as_S instead.", call. = FALSE)
}
}
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
# -- date
if (is.null(col_date)) {
2019-05-23 16:58:59 +02:00
col_date <- search_type_in_df(x = x, type = "date")
stop_if(is.null(col_date), "`col_date` must be set")
2019-01-15 12:45:24 +01:00
}
stop_ifnot(col_date %in% colnames(x),
"column `", col_date, "` not found")
2020-07-13 09:17:24 +02:00
# no grouped tibbles
2020-05-16 13:05:47 +02:00
x <- as.data.frame(x, stringsAsFactors = FALSE)
2018-08-10 15:01:05 +02:00
year <- function(x) {
2019-11-11 10:46:39 +01:00
# don't depend on lubridate or so, would be overkill for only this function
2019-10-11 17:21:02 +02:00
if (all(grepl("^[0-9]{4}$", x))) {
2020-06-17 21:35:10 +02:00
as.integer(x)
2018-08-10 15:01:05 +02:00
} else {
2019-10-11 17:21:02 +02:00
as.integer(format(as.Date(x), "%Y"))
2018-08-10 15:01:05 +02:00
}
}
2020-05-16 13:05:47 +02:00
df <- x
df[, col_ab] <- droplevels(as.rsi(df[, col_ab, drop = TRUE]))
if (I_as_S == TRUE) {
2020-05-16 13:05:47 +02:00
# then I as S
df[, col_ab] <- gsub("I", "S", df[, col_ab, drop = TRUE])
} else {
# then I as R
2020-05-16 13:05:47 +02:00
df[, col_ab] <- gsub("I", "R", df[, col_ab, drop = TRUE])
2018-08-10 15:01:05 +02:00
}
2020-05-16 13:05:47 +02:00
df[, col_ab] <- ifelse(is.na(df[, col_ab, drop = TRUE]), 0, df[, col_ab, drop = TRUE])
# remove rows with NAs
df <- subset(df, !is.na(df[, col_ab, drop = TRUE]))
df$year <- year(df[, col_date, drop = TRUE])
df <- as.data.frame(rbind(table(df[, c("year", col_ab)])), stringsAsFactors = FALSE)
df$year <- as.integer(rownames(df))
rownames(df) <- NULL
2020-07-13 09:17:24 +02:00
2020-05-16 13:05:47 +02:00
df <- subset(df, sum(df$R + df$S, na.rm = TRUE) >= minimum)
df_matrix <- as.matrix(df[, c("R", "S"), drop = FALSE])
stop_if(NROW(df) == 0, "there are no observations")
2018-08-10 15:01:05 +02:00
year_lowest <- min(df$year)
if (is.null(year_min)) {
year_min <- year_lowest
} else {
year_min <- max(year_min, year_lowest, na.rm = TRUE)
}
if (is.null(year_max)) {
2019-01-12 19:31:30 +01:00
year_max <- year(Sys.Date()) + 10
2018-08-10 15:01:05 +02:00
}
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
years <- list(year = seq(from = year_min, to = year_max, by = year_every))
2020-07-13 09:17:24 +02:00
2019-10-11 17:21:02 +02:00
if (model %in% c("binomial", "binom", "logit")) {
2019-01-15 12:45:24 +01:00
model <- "binomial"
model_lm <- with(df, glm(df_matrix ~ year, family = binomial))
2018-08-10 15:01:05 +02:00
if (info == TRUE) {
2019-10-11 17:21:02 +02:00
cat("\nLogistic regression model (logit) with binomial distribution")
cat("\n------------------------------------------------------------\n")
2019-01-15 12:45:24 +01:00
print(summary(model_lm))
2018-08-10 15:01:05 +02:00
}
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
predictmodel <- predict(model_lm, newdata = years, type = "response", se.fit = TRUE)
2018-08-10 15:01:05 +02:00
prediction <- predictmodel$fit
se <- predictmodel$se.fit
2020-07-13 09:17:24 +02:00
2019-10-11 17:21:02 +02:00
} else if (model %in% c("loglin", "poisson")) {
2019-01-15 12:45:24 +01:00
model <- "poisson"
model_lm <- with(df, glm(R ~ year, family = poisson))
2018-08-10 15:01:05 +02:00
if (info == TRUE) {
2019-10-11 17:21:02 +02:00
cat("\nLog-linear regression model (loglin) with poisson distribution")
cat("\n--------------------------------------------------------------\n")
2019-01-15 12:45:24 +01:00
print(summary(model_lm))
2018-08-10 15:01:05 +02:00
}
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
predictmodel <- predict(model_lm, newdata = years, type = "response", se.fit = TRUE)
2018-08-10 15:01:05 +02:00
prediction <- predictmodel$fit
se <- predictmodel$se.fit
2020-07-13 09:17:24 +02:00
2019-10-11 17:21:02 +02:00
} else if (model %in% c("lin", "linear")) {
2019-01-15 12:45:24 +01:00
model <- "linear"
model_lm <- with(df, lm((R / (R + S)) ~ year))
2018-08-10 15:01:05 +02:00
if (info == TRUE) {
2019-10-11 17:21:02 +02:00
cat("\nLinear regression model")
cat("\n-----------------------\n")
2019-01-15 12:45:24 +01:00
print(summary(model_lm))
2018-08-10 15:01:05 +02:00
}
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
predictmodel <- predict(model_lm, newdata = years, se.fit = TRUE)
2018-08-10 15:01:05 +02:00
prediction <- predictmodel$fit
se <- predictmodel$se.fit
2020-07-13 09:17:24 +02:00
2018-08-10 15:01:05 +02:00
} else {
stop("no valid model selected. See ?resistance_predict.")
2018-08-10 15:01:05 +02:00
}
2020-07-13 09:17:24 +02:00
2018-08-10 15:01:05 +02:00
# prepare the output dataframe
2019-01-15 12:45:24 +01:00
df_prediction <- data.frame(year = unlist(years),
value = prediction,
2020-05-16 13:05:47 +02:00
se_min = prediction - se,
se_max = prediction + se,
stringsAsFactors = FALSE)
2020-07-13 09:17:24 +02:00
2019-10-11 17:21:02 +02:00
if (model == "poisson") {
2020-05-16 13:05:47 +02:00
df_prediction$value <- as.integer(format(df_prediction$value, scientific = FALSE))
df_prediction$se_min <- as.integer(df_prediction$se_min)
df_prediction$se_max <- as.integer(df_prediction$se_max)
2018-08-10 15:01:05 +02:00
} else {
2020-05-16 13:05:47 +02:00
# se_max not above 1
df_prediction$se_max <- ifelse(df_prediction$se_max > 1, 1, df_prediction$se_max)
2018-08-10 15:01:05 +02:00
}
2020-05-16 13:05:47 +02:00
# se_min not below 0
df_prediction$se_min <- ifelse(df_prediction$se_min < 0, 0, df_prediction$se_min)
2020-07-13 09:17:24 +02:00
2020-05-16 13:05:47 +02:00
df_observations <- data.frame(year = df$year,
observations = df$R + df$S,
observed = df$R / (df$R + df$S),
stringsAsFactors = FALSE)
df_prediction <- df_prediction %pm>%
pm_left_join(df_observations, by = "year")
2020-05-16 13:05:47 +02:00
df_prediction$estimated <- df_prediction$value
2020-07-13 09:17:24 +02:00
2018-08-10 15:01:05 +02:00
if (preserve_measurements == TRUE) {
# replace estimated data by observed data
2020-05-16 13:05:47 +02:00
df_prediction$value <- ifelse(!is.na(df_prediction$observed), df_prediction$observed, df_prediction$value)
df_prediction$se_min <- ifelse(!is.na(df_prediction$observed), NA, df_prediction$se_min)
df_prediction$se_max <- ifelse(!is.na(df_prediction$observed), NA, df_prediction$se_max)
2018-08-10 15:01:05 +02:00
}
2020-07-13 09:17:24 +02:00
2020-05-16 13:05:47 +02:00
df_prediction$value <- ifelse(df_prediction$value > 1, 1, ifelse(df_prediction$value < 0, 0, df_prediction$value))
df_prediction <- df_prediction[order(df_prediction$year), ]
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
structure(
.Data = df_prediction,
class = c("resistance_predict", "data.frame"),
2019-05-13 16:35:48 +02:00
I_as_S = I_as_S,
2019-01-15 12:45:24 +01:00
model_title = model,
model = model_lm,
ab = col_ab
)
}
#' @rdname resistance_predict
#' @export
rsi_predict <- resistance_predict
2020-05-28 16:48:55 +02:00
#' @method plot resistance_predict
2019-01-15 12:45:24 +01:00
#' @export
#' @importFrom graphics plot axis arrows points
2019-01-15 12:45:24 +01:00
#' @rdname resistance_predict
2019-05-13 20:16:51 +02:00
plot.resistance_predict <- function(x, main = paste("Resistance Prediction of", x_name), ...) {
x_name <- paste0(ab_name(attributes(x)$ab), " (", attributes(x)$ab, ")")
meet_criteria(main, allow_class = "character", has_length = 1)
2020-07-13 09:17:24 +02:00
2019-05-13 16:35:48 +02:00
if (attributes(x)$I_as_S == TRUE) {
2019-01-15 12:45:24 +01:00
ylab <- "%R"
2019-05-13 16:35:48 +02:00
} else {
ylab <- "%IR"
2018-08-10 15:01:05 +02:00
}
2019-01-15 12:45:24 +01:00
plot(x = x$year,
y = x$value,
ylim = c(0, 1),
yaxt = "n", # no y labels
pch = 19, # closed dots
ylab = paste0("Percentage (", ylab, ")"),
xlab = "Year",
main = main,
2019-02-09 22:16:24 +01:00
sub = paste0("(n = ", sum(x$observations, na.rm = TRUE),
", model: ", attributes(x)$model_title, ")"),
cex.sub = 0.75)
2020-07-13 09:17:24 +02:00
2019-01-15 12:45:24 +01:00
axis(side = 2, at = seq(0, 1, 0.1), labels = paste0(0:10 * 10, "%"))
2020-07-13 09:17:24 +02:00
2019-02-11 10:27:10 +01:00
# hack for error bars: https://stackoverflow.com/a/22037078/4575331
2019-01-15 12:45:24 +01:00
arrows(x0 = x$year,
y0 = x$se_min,
x1 = x$year,
2019-02-11 10:27:10 +01:00
y1 = x$se_max,
length = 0.05, angle = 90, code = 3, lwd = 1.5)
2020-07-13 09:17:24 +02:00
2019-02-11 10:27:10 +01:00
# overlay grey points for prediction
2020-05-16 13:05:47 +02:00
points(x = subset(x, is.na(observations))$year,
y = subset(x, is.na(observations))$value,
2019-02-11 10:27:10 +01:00
pch = 19,
col = "grey40")
2018-08-10 15:01:05 +02:00
}
#' @rdname resistance_predict
#' @export
2019-02-11 10:27:10 +01:00
ggplot_rsi_predict <- function(x,
2019-05-13 20:16:51 +02:00
main = paste("Resistance Prediction of", x_name),
2019-02-11 10:27:10 +01:00
ribbon = TRUE,
...) {
x_name <- paste0(ab_name(attributes(x)$ab), " (", attributes(x)$ab, ")")
meet_criteria(main, allow_class = "character", has_length = 1)
meet_criteria(ribbon, allow_class = "logical", has_length = 1)
stop_ifnot_installed("ggplot2")
stop_ifnot(inherits(x, "resistance_predict"), "`x` must be a resistance prediction model created with resistance_predict()")
2020-07-13 09:17:24 +02:00
2019-05-13 16:35:48 +02:00
if (attributes(x)$I_as_S == TRUE) {
2019-01-15 12:45:24 +01:00
ylab <- "%R"
2019-05-13 16:35:48 +02:00
} else {
ylab <- "%IR"
2019-01-15 12:45:24 +01:00
}
2020-07-13 09:17:24 +02:00
2019-02-11 10:27:10 +01:00
p <- ggplot2::ggplot(x, ggplot2::aes(x = year, y = value)) +
2020-05-16 13:05:47 +02:00
ggplot2::geom_point(data = subset(x, !is.na(observations)),
2019-02-11 10:27:10 +01:00
size = 2) +
scale_y_percent(limits = c(0, 1)) +
ggplot2::labs(title = main,
y = paste0("Percentage (", ylab, ")"),
x = "Year",
caption = paste0("(n = ", sum(x$observations, na.rm = TRUE),
", model: ", attributes(x)$model_title, ")"))
2020-07-13 09:17:24 +02:00
2019-02-11 10:27:10 +01:00
if (ribbon == TRUE) {
p <- p + ggplot2::geom_ribbon(ggplot2::aes(ymin = se_min, ymax = se_max), alpha = 0.25)
} else {
p <- p + ggplot2::geom_errorbar(ggplot2::aes(ymin = se_min, ymax = se_max), na.rm = TRUE, width = 0.5)
}
p <- p +
# overlay grey points for prediction
2020-05-16 13:05:47 +02:00
ggplot2::geom_point(data = subset(x, is.na(observations)),
2019-02-11 10:27:10 +01:00
size = 2,
colour = "grey40")
p
2019-01-15 12:45:24 +01:00
}