1
0
mirror of https://github.com/msberends/AMR.git synced 2024-12-25 20:06:12 +01:00

resistance predict

This commit is contained in:
dr. M.S. (Matthijs) Berends 2019-01-15 12:45:24 +01:00
parent cda7087722
commit 046d195064
16 changed files with 306 additions and 209 deletions

View File

@ -1,6 +1,6 @@
Package: AMR Package: AMR
Version: 0.5.0.9009 Version: 0.5.0.9009
Date: 2019-01-12 Date: 2019-01-15
Title: Antimicrobial Resistance Analysis Title: Antimicrobial Resistance Analysis
Authors@R: c( Authors@R: c(
person( person(

View File

@ -20,6 +20,7 @@ S3method(kurtosis,default)
S3method(kurtosis,matrix) S3method(kurtosis,matrix)
S3method(plot,frequency_tbl) S3method(plot,frequency_tbl)
S3method(plot,mic) S3method(plot,mic)
S3method(plot,resistance_predict)
S3method(plot,rsi) S3method(plot,rsi)
S3method(print,atc) S3method(print,atc)
S3method(print,frequency_tbl) S3method(print,frequency_tbl)
@ -76,6 +77,7 @@ export(g.test)
export(geom_rsi) export(geom_rsi)
export(get_locale) export(get_locale)
export(ggplot_rsi) export(ggplot_rsi)
export(ggplot_rsi_predict)
export(guess_ab_col) export(guess_ab_col)
export(guess_atc) export(guess_atc)
export(guess_mo) export(guess_mo)
@ -215,6 +217,7 @@ importFrom(dplyr,mutate_all)
importFrom(dplyr,mutate_at) importFrom(dplyr,mutate_at)
importFrom(dplyr,n) importFrom(dplyr,n)
importFrom(dplyr,n_distinct) importFrom(dplyr,n_distinct)
importFrom(dplyr,n_groups)
importFrom(dplyr,progress_estimated) importFrom(dplyr,progress_estimated)
importFrom(dplyr,pull) importFrom(dplyr,pull)
importFrom(dplyr,row_number) importFrom(dplyr,row_number)
@ -228,6 +231,7 @@ importFrom(dplyr,top_n)
importFrom(dplyr,ungroup) importFrom(dplyr,ungroup)
importFrom(dplyr,vars) importFrom(dplyr,vars)
importFrom(grDevices,boxplot.stats) importFrom(grDevices,boxplot.stats)
importFrom(graphics,arrows)
importFrom(graphics,axis) importFrom(graphics,axis)
importFrom(graphics,barplot) importFrom(graphics,barplot)
importFrom(graphics,hist) importFrom(graphics,hist)

View File

@ -12,6 +12,12 @@
* Function `mo_renamed()` to get a list of all returned values from `as.mo()` that have had taxonomic renaming * Function `mo_renamed()` to get a list of all returned values from `as.mo()` that have had taxonomic renaming
* Function `age()` to calculate the (patients) age in years * Function `age()` to calculate the (patients) age in years
* Function `age_groups()` to split ages into custom or predefined groups (like children or elderly). This allows for easier demographic antimicrobial resistance analysis per age group. * Function `age_groups()` to split ages into custom or predefined groups (like children or elderly). This allows for easier demographic antimicrobial resistance analysis per age group.
* Function `ggplot_rsi_predict()` as well as the base R `plot()` function can now be used for resistance prediction calculated with `resistance_predict()`:
```r
x <- resistance_predict(septic_patients, col_ab = "amox")
plot(x)
ggplot_rsi_predict(x)
```
* Functions `filter_first_isolate()` and `filter_first_weighted_isolate()` to shorten and fasten filtering on data sets with antimicrobial results, e.g.: * Functions `filter_first_isolate()` and `filter_first_weighted_isolate()` to shorten and fasten filtering on data sets with antimicrobial results, e.g.:
```r ```r
septic_patients %>% filter_first_isolate(...) septic_patients %>% filter_first_isolate(...)

View File

@ -231,9 +231,8 @@ eucast_rules <- function(tbl,
# try to find columns based on type # try to find columns based on type
# -- mo # -- mo
if (is.null(col_mo) & "mo" %in% lapply(tbl, class)) { if (is.null(col_mo)) {
col_mo <- colnames(tbl)[lapply(tbl, class) == "mo"][1] col_mo <- search_type_in_df(tbl = tbl, type = "mo")
message(blue(paste0("NOTE: Using column `", bold(col_mo), "` as input for `col_mo`.")))
} }
if (is.null(col_mo)) { if (is.null(col_mo)) {
stop("`col_mo` must be set.", call. = FALSE) stop("`col_mo` must be set.", call. = FALSE)

View File

@ -23,7 +23,7 @@
#' #'
#' Determine first (weighted) isolates of all microorganisms of every patient per episode and (if needed) per specimen type. #' Determine first (weighted) isolates of all microorganisms of every patient per episode and (if needed) per specimen type.
#' @param tbl a \code{data.frame} containing isolates. #' @param tbl a \code{data.frame} containing isolates.
#' @param col_date column name of the result date (or date that is was received on the lab), defaults to the first column of class \code{Date} #' @param col_date column name of the result date (or date that is was received on the lab), defaults to the first column of with a date class
#' @param col_patient_id column name of the unique IDs of the patients, defaults to the first column that starts with 'patient' or 'patid' (case insensitive) #' @param col_patient_id column name of the unique IDs of the patients, defaults to the first column that starts with 'patient' or 'patid' (case insensitive)
#' @param col_mo column name of the unique IDs of the microorganisms (see \code{\link{mo}}), defaults to the first column of class \code{mo}. Values will be coerced using \code{\link{as.mo}}. #' @param col_mo column name of the unique IDs of the microorganisms (see \code{\link{mo}}), defaults to the first column of class \code{mo}. Values will be coerced using \code{\link{as.mo}}.
#' @param col_testcode column name of the test codes. Use \code{col_testcode = NULL} to \strong{not} exclude certain test codes (like test codes for screening). In that case \code{testcodes_exclude} will be ignored. #' @param col_testcode column name of the test codes. Use \code{col_testcode = NULL} to \strong{not} exclude certain test codes (like test codes for screening). In that case \code{testcodes_exclude} will be ignored.
@ -187,9 +187,8 @@ first_isolate <- function(tbl,
# try to find columns based on type # try to find columns based on type
# -- mo # -- mo
if (is.null(col_mo) & "mo" %in% lapply(tbl, class)) { if (is.null(col_mo)) {
col_mo <- colnames(tbl)[lapply(tbl, class) == "mo"][1] col_mo <- search_type_in_df(tbl = tbl, type = "mo")
message(blue(paste0("NOTE: Using column `", bold(col_mo), "` as input for `col_mo`.")))
} }
if (is.null(col_mo)) { if (is.null(col_mo)) {
stop("`col_mo` must be set.", call. = FALSE) stop("`col_mo` must be set.", call. = FALSE)
@ -197,33 +196,25 @@ first_isolate <- function(tbl,
# -- date # -- date
if (is.null(col_date)) { if (is.null(col_date)) {
for (i in 1:ncol(tbl)) { col_date <- search_type_in_df(tbl = tbl, type = "date")
if ("Date" %in% class(tbl %>% pull(i)) | "POSIXct" %in% class(tbl %>% pull(i))) {
col_date <- colnames(tbl)[i]
message(blue(paste0("NOTE: Using column `", bold(col_date), "` as input for `col_date`.")))
break
}
}
} }
if (is.null(col_date)) { if (is.null(col_date)) {
stop("`col_date` must be set.", call. = FALSE) stop("`col_date` must be set.", call. = FALSE)
} }
# convert to Date (pipes for supporting tibbles too) # convert to Date (pipes/pull for supporting tibbles too)
tbl[, col_date] <- tbl %>% pull(col_date) %>% as.Date() tbl[, col_date] <- tbl %>% pull(col_date) %>% as.Date()
# -- patient id # -- patient id
if (is.null(col_patient_id) & any(colnames(tbl) %like% "^(patient|patid)")) { if (is.null(col_patient_id)) {
col_patient_id <- colnames(tbl)[colnames(tbl) %like% "^(patient|patid)"][1] col_patient_id <- search_type_in_df(tbl = tbl, type = "patient_id")
message(blue(paste0("NOTE: Using column `", bold(col_patient_id), "` as input for `col_patient_id`.")))
} }
if (is.null(col_patient_id)) { if (is.null(col_patient_id)) {
stop("`col_patient_id` must be set.", call. = FALSE) stop("`col_patient_id` must be set.", call. = FALSE)
} }
# -- key antibiotics # -- key antibiotics
if (is.null(col_keyantibiotics) & any(colnames(tbl) %like% "^key.*(ab|antibiotics)")) { if (is.null(col_keyantibiotics)) {
col_keyantibiotics <- colnames(tbl)[colnames(tbl) %like% "^key.*(ab|antibiotics)"][1] col_keyantibiotics <- search_type_in_df(tbl = tbl, type = "keyantibiotics")
message(blue(paste0("NOTE: Using column `", bold(col_keyantibiotics), "` as input for `col_keyantibiotics`. Use ", bold("col_keyantibiotics = FALSE"), " to prevent this.")))
} }
if (isFALSE(col_keyantibiotics)) { if (isFALSE(col_keyantibiotics)) {
col_keyantibiotics <- NULL col_keyantibiotics <- NULL

View File

@ -101,9 +101,8 @@ key_antibiotics <- function(tbl,
# try to find columns based on type # try to find columns based on type
# -- mo # -- mo
if (is.null(col_mo) & "mo" %in% lapply(tbl, class)) { if (is.null(col_mo)) {
col_mo <- colnames(tbl)[lapply(tbl, class) == "mo"][1] col_mo <- search_type_in_df(tbl = tbl, type = "mo")
message(blue(paste0("NOTE: Using column `", bold(col_mo), "` as input for `col_mo`.")))
} }
if (is.null(col_mo)) { if (is.null(col_mo)) {
stop("`col_mo` must be set.", call. = FALSE) stop("`col_mo` must be set.", call. = FALSE)
@ -114,7 +113,6 @@ key_antibiotics <- function(tbl,
GramPos_1, GramPos_2, GramPos_3, GramPos_4, GramPos_5, GramPos_6, GramPos_1, GramPos_2, GramPos_3, GramPos_4, GramPos_5, GramPos_6,
GramNeg_1, GramNeg_2, GramNeg_3, GramNeg_4, GramNeg_5, GramNeg_6) GramNeg_1, GramNeg_2, GramNeg_3, GramNeg_4, GramNeg_5, GramNeg_6)
col.list <- check_available_columns(tbl = tbl, col.list = col.list, info = warnings) col.list <- check_available_columns(tbl = tbl, col.list = col.list, info = warnings)
print(col.list)
universal_1 <- col.list[universal_1] universal_1 <- col.list[universal_1]
universal_2 <- col.list[universal_2] universal_2 <- col.list[universal_2]
universal_3 <- col.list[universal_3] universal_3 <- col.list[universal_3]

View File

@ -113,9 +113,8 @@ mdro <- function(tbl,
# try to find columns based on type # try to find columns based on type
# -- mo # -- mo
if (is.null(col_mo) & "mo" %in% lapply(tbl, class)) { if (is.null(col_mo)) {
col_mo <- colnames(tbl)[lapply(tbl, class) == "mo"][1] col_mo <- search_type_in_df(tbl = tbl, type = "mo")
message(blue(paste0("NOTE: Using column `", bold(col_mo), "` as input for `col_mo`.")))
} }
if (is.null(col_mo)) { if (is.null(col_mo)) {
stop("`col_mo` must be set.", call. = FALSE) stop("`col_mo` must be set.", call. = FALSE)

View File

@ -123,3 +123,48 @@ size_humanreadable <- function(bytes, decimals = 1) {
out <- paste(sprintf(paste0("%.", decimals, "f"), bytes / (1024 ^ factor)), size[factor + 1]) out <- paste(sprintf(paste0("%.", decimals, "f"), bytes / (1024 ^ factor)), size[factor + 1])
out out
} }
#' @importFrom crayon blue bold
#' @importFrom dplyr %>% pull
search_type_in_df <- function(tbl, type) {
# try to find columns based on type
found <- NULL
# -- mo
if (type == "mo") {
if ("mo" %in% lapply(tbl, class)) {
found <- colnames(tbl)[lapply(tbl, class) == "mo"][1]
}
}
# -- key antibiotics
if (type == "keyantibiotics") {
if (any(colnames(tbl) %like% "^key.*(ab|antibiotics)")) {
found <- colnames(tbl)[colnames(tbl) %like% "^key.*(ab|antibiotics)"][1]
}
}
# -- date
if (type == "date") {
for (i in 1:ncol(tbl)) {
if ("Date" %in% class(tbl %>% pull(i)) | "POSIXct" %in% class(tbl %>% pull(i))) {
found <- colnames(tbl)[i]
break
}
}
}
# -- patient id
if (type == "patient_id") {
if (any(colnames(tbl) %like% "^(patient|patid)")) {
found <- colnames(tbl)[colnames(tbl) %like% "^(patient|patid)"][1]
}
}
if (!is.null(found)) {
msg <- paste0("NOTE: Using column `", bold(found), "` as input for `col_", type, "`.")
if (type == "keyantibiotics") {
msg <- paste(msg, "Use", bold("col_keyantibiotics = FALSE"), "to prevent this.")
}
message(blue(msg))
}
found
}

View File

@ -23,17 +23,18 @@
#' #'
#' Create a prediction model to predict antimicrobial resistance for the next years on statistical solid ground. Standard errors (SE) will be returned as columns \code{se_min} and \code{se_max}. See Examples for a real live example. #' Create a prediction model to predict antimicrobial resistance for the next years on statistical solid ground. Standard errors (SE) will be returned as columns \code{se_min} and \code{se_max}. See Examples for a real live example.
#' @inheritParams first_isolate #' @inheritParams first_isolate
#' @inheritParams graphics::plot
#' @param col_ab column name of \code{tbl} with antimicrobial interpretations (\code{R}, \code{I} and \code{S}) #' @param col_ab column name of \code{tbl} with antimicrobial interpretations (\code{R}, \code{I} and \code{S})
#' @param col_date column name of the date, will be used to calculate years if this column doesn't consist of years already #' @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 \code{col_date} #' @param year_min lowest year to use in the prediction model, dafaults to the lowest year in \code{col_date}
#' @param year_max highest year to use in the prediction model, defaults to 10 years after today #' @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 \code{year_max} #' @param year_every unit of sequence between lowest year found in the data and \code{year_max}
#' @param minimum minimal amount of available isolates per year to include. Years containing less observations will be estimated by the model. #' @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. Valid values are \code{"binomial"} (or \code{"binom"} or \code{"logit"}) or \code{"loglin"} or \code{"linear"} (or \code{"lin"}). #' @param model the statistical model of choice. Valid values are \code{"binomial"} (or \code{"binom"} or \code{"logit"}) or \code{"loglin"} (or \code{"poisson"}) or \code{"linear"} (or \code{"lin"}).
#' @param I_as_R a logical to indicate whether values \code{I} should be treated as \code{R} #' @param I_as_R a logical to indicate whether values \code{I} should be treated as \code{R}
#' @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 \code{NA}. #' @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 \code{NA}.
#' @param info a logical to indicate whether textual analysis should be printed with the name and \code{\link{summary}} of the statistical model. #' @param info a logical to indicate whether textual analysis should be printed with the name and \code{\link{summary}} of the statistical model.
#' @return \code{data.frame} with columns: #' @return \code{data.frame} with extra class \code{"resistance_predict"} with columns:
#' \itemize{ #' \itemize{
#' \item{\code{year}} #' \item{\code{year}}
#' \item{\code{value}, the same as \code{estimated} when \code{preserve_measurements = FALSE}, and a combination of \code{observed} and \code{estimated} otherwise} #' \item{\code{value}, the same as \code{estimated} when \code{preserve_measurements = FALSE}, and a combination of \code{observed} and \code{estimated} otherwise}
@ -47,42 +48,23 @@
#' @rdname resistance_predict #' @rdname resistance_predict
#' @export #' @export
#' @importFrom stats predict glm lm #' @importFrom stats predict glm lm
#' @importFrom dplyr %>% pull mutate mutate_at n group_by_at summarise filter filter_at all_vars n_distinct arrange case_when #' @importFrom dplyr %>% pull mutate mutate_at n group_by_at summarise filter filter_at all_vars n_distinct arrange case_when n_groups
#' @inheritSection AMR Read more on our website! #' @inheritSection AMR Read more on our website!
#' @examples #' @examples
#' \dontrun{ #' x <- resistance_predict(septic_patients, col_ab = "amox", year_min = 2010)
#' # use it with base R: #' plot(x)
#' resistance_predict(tbl = tbl[which(first_isolate == TRUE & genus == "Haemophilus"),], #' ggplot_rsi_predict(x)
#' col_ab = "amcl", col_date = "date")
#' #'
#' # or use dplyr so you can actually read it: #' # use dplyr so you can actually read it:
#' library(dplyr) #' library(dplyr)
#' tbl %>% #' x <- septic_patients %>%
#' filter(first_isolate == TRUE, #' filter_first_isolate() %>%
#' genus == "Haemophilus") %>% #' filter(mo_genus(mo) == "Staphylococcus") %>%
#' resistance_predict(amcl, date) #' resistance_predict("peni")
#' } #' plot(x)
#' #'
#' #'
#' # real live example: #' # create nice plots with ggplot yourself
#' library(dplyr)
#' septic_patients %>%
#' # get bacteria properties like genus and species
#' left_join_microorganisms("mo") %>%
#' # calculate first isolates
#' mutate(first_isolate = first_isolate(.)) %>%
#' # filter on first E. coli isolates
#' filter(genus == "Escherichia",
#' species == "coli",
#' first_isolate == TRUE) %>%
#' # predict resistance of cefotaxime for next years
#' resistance_predict(col_ab = "cfot",
#' col_date = "date",
#' year_max = 2025,
#' preserve_measurements = TRUE,
#' minimum = 0)
#'
#' # create nice plots with ggplot
#' if (!require(ggplot2)) { #' if (!require(ggplot2)) {
#' #'
#' data <- septic_patients %>% #' data <- septic_patients %>%
@ -90,7 +72,7 @@
#' resistance_predict(col_ab = "amox", #' resistance_predict(col_ab = "amox",
#' col_date = "date", #' col_date = "date",
#' info = FALSE, #' info = FALSE,
#' minimum = 15) #' minimum = 15)
#' #'
#' ggplot(data, #' ggplot(data,
#' aes(x = year)) + #' aes(x = year)) +
@ -110,7 +92,7 @@
#' } #' }
resistance_predict <- function(tbl, resistance_predict <- function(tbl,
col_ab, col_ab,
col_date, col_date = NULL,
year_min = NULL, year_min = NULL,
year_max = NULL, year_max = NULL,
year_every = 1, year_every = 1,
@ -128,23 +110,23 @@ resistance_predict <- function(tbl,
stop('Column ', col_ab, ' not found.') stop('Column ', col_ab, ' not found.')
} }
# -- date
if (is.null(col_date)) {
col_date <- search_type_in_df(tbl = tbl, type = "date")
}
if (is.null(col_date)) {
stop("`col_date` must be set.", call. = FALSE)
}
if (!col_date %in% colnames(tbl)) { if (!col_date %in% colnames(tbl)) {
stop('Column ', col_date, ' not found.') stop('Column ', col_date, ' not found.')
} }
if ('grouped_df' %in% class(tbl)) {
if (n_groups(tbl) > 1) {
# no grouped tibbles please, mutate will throw errors # no grouped tibbles please, mutate will throw errors
tbl <- base::as.data.frame(tbl, stringsAsFactors = FALSE) tbl <- base::as.data.frame(tbl, stringsAsFactors = FALSE)
} }
if (I_as_R == TRUE) {
tbl[, col_ab] <- gsub('I', 'R', tbl %>% pull(col_ab))
}
tbl <- tbl %>%
mutate_at(col_ab, as.rsi) %>%
filter_at(col_ab, all_vars(!is.na(.)))
tbl[, col_ab] <- droplevels(tbl[, col_ab])
year <- function(x) { year <- function(x) {
if (all(grepl('^[0-9]{4}$', x))) { if (all(grepl('^[0-9]{4}$', x))) {
x x
@ -154,13 +136,23 @@ resistance_predict <- function(tbl,
} }
df <- tbl %>% df <- tbl %>%
mutate(year = tbl %>% pull(col_date) %>% year()) %>% mutate_at(col_ab, as.rsi) %>%
mutate_at(col_ab, droplevels) %>%
mutate_at(col_ab, funs(
if (I_as_R == TRUE) {
gsub("I", "R", .)
} else {
gsub("I", "S", .)
}
)) %>%
filter_at(col_ab, all_vars(!is.na(.))) %>%
mutate(year = pull(., col_date) %>% year()) %>%
group_by_at(c('year', col_ab)) %>% group_by_at(c('year', col_ab)) %>%
summarise(n()) summarise(n())
if (df %>% pull(col_ab) %>% n_distinct(na.rm = TRUE) < 2) { if (df %>% pull(col_ab) %>% n_distinct(na.rm = TRUE) < 2) {
stop("No variety in antimicrobial interpretations - all isolates are '", stop("No variety in antimicrobial interpretations - all isolates are '",
df %>% pull(col_ab) %>% unique() %>% .[!is.na(.)], "'.", df %>% pull(col_ab) %>% unique(), "'.",
call. = FALSE) call. = FALSE)
} }
@ -168,8 +160,11 @@ resistance_predict <- function(tbl,
df <- df %>% df <- df %>%
filter(!is.na(antibiotic)) %>% filter(!is.na(antibiotic)) %>%
tidyr::spread(antibiotic, observations, fill = 0) %>% tidyr::spread(antibiotic, observations, fill = 0) %>%
mutate(total = R + S) %>% filter((R + S) >= minimum)
filter(total >= minimum) df_matrix <- df %>%
ungroup() %>%
select(R, S) %>%
as.matrix()
if (NROW(df) == 0) { if (NROW(df) == 0) {
stop('There are no observations.') stop('There are no observations.')
@ -185,41 +180,44 @@ resistance_predict <- function(tbl,
year_max <- year(Sys.Date()) + 10 year_max <- year(Sys.Date()) + 10
} }
years_predict <- seq(from = year_min, to = year_max, by = year_every) years <- list(year = seq(from = year_min, to = year_max, by = year_every))
if (model %in% c('binomial', 'binom', 'logit')) { if (model %in% c('binomial', 'binom', 'logit')) {
logitmodel <- with(df, glm(cbind(R, S) ~ year, family = binomial)) model <- "binomial"
model_lm <- with(df, glm(df_matrix ~ year, family = binomial))
if (info == TRUE) { if (info == TRUE) {
cat('\nLogistic regression model (logit) with binomial distribution') cat('\nLogistic regression model (logit) with binomial distribution')
cat('\n------------------------------------------------------------\n') cat('\n------------------------------------------------------------\n')
print(summary(logitmodel)) print(summary(model_lm))
} }
predictmodel <- predict(logitmodel, newdata = with(df, list(year = years_predict)), type = "response", se.fit = TRUE) predictmodel <- predict(model_lm, newdata = years, type = "response", se.fit = TRUE)
prediction <- predictmodel$fit prediction <- predictmodel$fit
se <- predictmodel$se.fit se <- predictmodel$se.fit
} else if (model == 'loglin') { } else if (model %in% c('loglin', 'poisson')) {
loglinmodel <- with(df, glm(R ~ year, family = poisson)) model <- "poisson"
model_lm <- with(df, glm(R ~ year, family = poisson))
if (info == TRUE) { if (info == TRUE) {
cat('\nLog-linear regression model (loglin) with poisson distribution') cat('\nLog-linear regression model (loglin) with poisson distribution')
cat('\n--------------------------------------------------------------\n') cat('\n--------------------------------------------------------------\n')
print(summary(loglinmodel)) print(summary(model_lm))
} }
predictmodel <- predict(loglinmodel, newdata = with(df, list(year = years_predict)), type = "response", se.fit = TRUE) predictmodel <- predict(model_lm, newdata = years, type = "response", se.fit = TRUE)
prediction <- predictmodel$fit prediction <- predictmodel$fit
se <- predictmodel$se.fit se <- predictmodel$se.fit
} else if (model %in% c('lin', 'linear')) { } else if (model %in% c('lin', 'linear')) {
linmodel <- with(df, lm((R / (R + S)) ~ year)) model <- "linear"
model_lm <- with(df, lm((R / (R + S)) ~ year))
if (info == TRUE) { if (info == TRUE) {
cat('\nLinear regression model') cat('\nLinear regression model')
cat('\n-----------------------\n') cat('\n-----------------------\n')
print(summary(linmodel)) print(summary(model_lm))
} }
predictmodel <- predict(linmodel, newdata = with(df, list(year = years_predict)), se.fit = TRUE) predictmodel <- predict(model_lm, newdata = years, se.fit = TRUE)
prediction <- predictmodel$fit prediction <- predictmodel$fit
se <- predictmodel$se.fit se <- predictmodel$se.fit
@ -228,64 +226,117 @@ resistance_predict <- function(tbl,
} }
# prepare the output dataframe # prepare the output dataframe
prediction <- data.frame(year = years_predict, value = prediction, stringsAsFactors = FALSE) df_prediction <- data.frame(year = unlist(years),
value = prediction,
stringsAsFactors = FALSE) %>%
prediction$se_min <- prediction$value - se mutate(se_min = value - se,
prediction$se_max <- prediction$value + se se_max = value + se)
if (model == 'loglin') { if (model == 'poisson') {
prediction$value <- prediction$value %>% df_prediction <- df_prediction %>%
format(scientific = FALSE) %>% mutate(value = value %>%
as.integer() format(scientific = FALSE) %>%
prediction$se_min <- prediction$se_min %>% as.integer() as.integer(),
prediction$se_max <- prediction$se_max %>% as.integer() se_min = as.integer(se_min),
se_max = as.integer(se_max))
colnames(prediction) <- c('year', 'amountR', 'se_max', 'se_min')
} else { } else {
prediction$se_max[which(prediction$se_max > 1)] <- 1 df_prediction <- df_prediction %>%
# se_max not above 1
mutate(se_max = ifelse(se_max > 1, 1, se_max))
} }
prediction$se_min[which(prediction$se_min < 0)] <- 0 df_prediction <- df_prediction %>%
prediction$observations = NA # se_min not below 0
mutate(se_min = ifelse(se_min < 0, 0, se_min))
total <- prediction df_observations <- df %>%
ungroup() %>%
transmute(year,
observations = R + S,
observed = R / (R + S))
df_prediction <- df_prediction %>%
left_join(df_observations, by = "year") %>%
mutate(estimated = value)
if (preserve_measurements == TRUE) { if (preserve_measurements == TRUE) {
# replace estimated data by observed data # replace estimated data by observed data
if (I_as_R == TRUE) { df_prediction <- df_prediction %>%
if (!'I' %in% colnames(df)) { mutate(value = ifelse(!is.na(observed), observed, value),
df$I <- 0 se_min = ifelse(!is.na(observed), NA, se_min),
} se_max = ifelse(!is.na(observed), NA, se_max))
df$value <- df$R / rowSums(df[, c('R', 'S', 'I')])
} else {
df$value <- df$R / rowSums(df[, c('R', 'S')])
}
measurements <- data.frame(year = df$year,
value = df$value,
se_min = NA,
se_max = NA,
observations = df$total,
stringsAsFactors = FALSE)
colnames(measurements) <- colnames(prediction)
total <- rbind(measurements,
prediction %>% filter(!year %in% df$year))
if (model %in% c('binomial', 'binom', 'logit')) {
total <- total %>% mutate(observed = ifelse(is.na(observations), NA, value),
estimated = prediction$value)
}
} }
if ("value" %in% colnames(total)) { df_prediction <- df_prediction %>%
total <- total %>% mutate(value = case_when(value > 1 ~ 1,
mutate(value = case_when(value > 1 ~ 1, value < 0 ~ 0,
value < 0 ~ 0, TRUE ~ value)) %>%
TRUE ~ value)) arrange(year)
}
total %>% arrange(year)
structure(
.Data = df_prediction,
class = c("resistance_predict", "data.frame"),
I_as_R = I_as_R,
model_title = model,
model = model_lm,
ab = col_ab
)
} }
#' @rdname resistance_predict #' @rdname resistance_predict
#' @export #' @export
rsi_predict <- resistance_predict rsi_predict <- resistance_predict
#' @exportMethod plot.mic
#' @export
#' @importFrom dplyr %>% group_by summarise
#' @importFrom graphics plot axis arrows
#' @rdname resistance_predict
plot.resistance_predict <- function(x, main = paste("Resistance prediction of", attributes(x)$ab), ...) {
if (attributes(x)$I_as_R == TRUE) {
ylab <- "%IR"
} else {
ylab <- "%R"
}
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,
sub = paste0("(model: ", attributes(x)$model_title, ")"))
axis(side = 2, at = seq(0, 1, 0.1), labels = paste0(0:10 * 10, "%"))
# arrows hack: https://stackoverflow.com/a/22037078/4575331
arrows(x0 = x$year,
y0 = x$se_min,
x1 = x$year,
y1 = x$se_max, length = 0.05, angle = 90, code = 3)
}
#' @rdname resistance_predict
#' @export
ggplot_rsi_predict <- function(x, main = paste("Resistance prediction of", attributes(x)$ab), ...) {
if (!"resistance_predict" %in% class(x)) {
stop("`x` must be a resistance prediction model created with resistance_predict().")
}
if (attributes(x)$I_as_R == TRUE) {
ylab <- "%IR"
} else {
ylab <- "%R"
}
suppressWarnings(
ggplot2::ggplot(x, ggplot2::aes(x = year, y = value)) +
ggplot2::geom_col() +
ggplot2::geom_errorbar(ggplot2::aes(ymin = se_min, ymax = se_max)) +
scale_y_percent() +
labs(title = main,
y = paste0("Percentage (", ylab, ")"),
x = "Year",
caption = paste0("(model: ", attributes(x)$model_title, ")"))
)
}

View File

@ -238,6 +238,8 @@
<li>Function <code><a href="../reference/mo_renamed.html">mo_renamed()</a></code> to get a list of all returned values from <code><a href="../reference/as.mo.html">as.mo()</a></code> that have had taxonomic renaming</li> <li>Function <code><a href="../reference/mo_renamed.html">mo_renamed()</a></code> to get a list of all returned values from <code><a href="../reference/as.mo.html">as.mo()</a></code> that have had taxonomic renaming</li>
<li>Function <code><a href="../reference/age.html">age()</a></code> to calculate the (patients) age in years</li> <li>Function <code><a href="../reference/age.html">age()</a></code> to calculate the (patients) age in years</li>
<li>Function <code><a href="../reference/age_groups.html">age_groups()</a></code> to split ages into custom or predefined groups (like children or elderly). This allows for easier demographic antimicrobial resistance analysis per age group.</li> <li>Function <code><a href="../reference/age_groups.html">age_groups()</a></code> to split ages into custom or predefined groups (like children or elderly). This allows for easier demographic antimicrobial resistance analysis per age group.</li>
<li>Function <code><a href="../reference/resistance_predict.html">ggplot_rsi_predict()</a></code> as well as the base R <code><a href="https://www.rdocumentation.org/packages/graphics/topics/plot">plot()</a></code> function can now be used for resistance prediction calculated with <code><a href="../reference/resistance_predict.html">resistance_predict()</a></code>: <code>r x &lt;- resistance_predict(septic_patients, col_ab = "amox") plot(x) ggplot_rsi_predict(x)</code>
</li>
<li>Functions <code><a href="../reference/first_isolate.html">filter_first_isolate()</a></code> and <code><a href="../reference/first_isolate.html">filter_first_weighted_isolate()</a></code> to shorten and fasten filtering on data sets with antimicrobial results, e.g.: <code>r septic_patients %&gt;% filter_first_isolate(...) # or filter_first_isolate(septic_patients, ...)</code> is equal to: <code>r septic_patients %&gt;% mutate(only_firsts = first_isolate(septic_patients, ...)) %&gt;% filter(only_firsts == TRUE) %&gt;% select(-only_firsts)</code> <li>Functions <code><a href="../reference/first_isolate.html">filter_first_isolate()</a></code> and <code><a href="../reference/first_isolate.html">filter_first_weighted_isolate()</a></code> to shorten and fasten filtering on data sets with antimicrobial results, e.g.: <code>r septic_patients %&gt;% filter_first_isolate(...) # or filter_first_isolate(septic_patients, ...)</code> is equal to: <code>r septic_patients %&gt;% mutate(only_firsts = first_isolate(septic_patients, ...)) %&gt;% filter(only_firsts == TRUE) %&gt;% select(-only_firsts)</code>
</li> </li>
<li>New vignettes about how to conduct AMR analysis, predict antimicrobial resistance, use the <em>G</em>-test and more. These are also available (and even easier readable) on our website: <a href="https://msberends.gitlab.io/AMR" class="uri">https://msberends.gitlab.io/AMR</a>.</li> <li>New vignettes about how to conduct AMR analysis, predict antimicrobial resistance, use the <em>G</em>-test and more. These are also available (and even easier readable) on our website: <a href="https://msberends.gitlab.io/AMR" class="uri">https://msberends.gitlab.io/AMR</a>.</li>

View File

@ -250,7 +250,7 @@
</tr> </tr>
<tr> <tr>
<th>col_date</th> <th>col_date</th>
<td><p>column name of the result date (or date that is was received on the lab), defaults to the first column of class <code>Date</code></p></td> <td><p>column name of the result date (or date that is was received on the lab), defaults to the first column of with a date class</p></td>
</tr> </tr>
<tr> <tr>
<th>col_patient_id</th> <th>col_patient_id</th>

View File

@ -415,7 +415,7 @@
</tr><tr> </tr><tr>
<td> <td>
<p><code><a href="resistance_predict.html">resistance_predict()</a></code> <code><a href="resistance_predict.html">rsi_predict()</a></code> </p> <p><code><a href="resistance_predict.html">resistance_predict()</a></code> <code><a href="resistance_predict.html">rsi_predict()</a></code> <code><a href="resistance_predict.html">plot(<i>&lt;resistance_predict&gt;</i>)</a></code> <code><a href="resistance_predict.html">ggplot_rsi_predict()</a></code> </p>
</td> </td>
<td><p>Predict antimicrobial resistance</p></td> <td><p>Predict antimicrobial resistance</p></td>
</tr><tr> </tr><tr>

View File

@ -227,14 +227,22 @@
</div> </div>
<pre class="usage"><span class='fu'>resistance_predict</span>(<span class='no'>tbl</span>, <span class='no'>col_ab</span>, <span class='no'>col_date</span>, <span class='kw'>year_min</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <pre class="usage"><span class='fu'>resistance_predict</span>(<span class='no'>tbl</span>, <span class='no'>col_ab</span>, <span class='kw'>col_date</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>year_min</span> <span class='kw'>=</span> <span class='kw'>NULL</span>,
<span class='kw'>year_max</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>year_every</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>30</span>, <span class='kw'>year_max</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>year_every</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>30</span>,
<span class='kw'>model</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>, <span class='kw'>I_as_R</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>preserve_measurements</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>model</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>, <span class='kw'>I_as_R</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>preserve_measurements</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>,
<span class='kw'>info</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) <span class='kw'>info</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
<span class='fu'>rsi_predict</span>(<span class='no'>tbl</span>, <span class='no'>col_ab</span>, <span class='no'>col_date</span>, <span class='kw'>year_min</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>year_max</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='fu'>rsi_predict</span>(<span class='no'>tbl</span>, <span class='no'>col_ab</span>, <span class='kw'>col_date</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>year_min</span> <span class='kw'>=</span> <span class='kw'>NULL</span>,
<span class='kw'>year_every</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>30</span>, <span class='kw'>model</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>, <span class='kw'>I_as_R</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>year_max</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>year_every</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>30</span>,
<span class='kw'>preserve_measurements</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>info</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)</pre> <span class='kw'>model</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>, <span class='kw'>I_as_R</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>preserve_measurements</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>,
<span class='kw'>info</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
<span class='co'># S3 method for resistance_predict</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/graphics/topics/plot'>plot</a></span>(<span class='no'>x</span>,
<span class='kw'>main</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/paste'>paste</a></span>(<span class='st'>"Resistance prediction of"</span>, <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/attributes'>attributes</a></span>(<span class='no'>x</span>)$<span class='no'>ab</span>), <span class='no'>...</span>)
<span class='fu'>ggplot_rsi_predict</span>(<span class='no'>x</span>, <span class='kw'>main</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/paste'>paste</a></span>(<span class='st'>"Resistance prediction of"</span>,
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/attributes'>attributes</a></span>(<span class='no'>x</span>)$<span class='no'>ab</span>), <span class='no'>...</span>)</pre>
<h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a>Arguments</h2> <h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a>Arguments</h2>
<table class="ref-arguments"> <table class="ref-arguments">
@ -249,7 +257,7 @@
</tr> </tr>
<tr> <tr>
<th>col_date</th> <th>col_date</th>
<td><p>column name of the date, will be used to calculate years if this column doesn't consist of years already</p></td> <td><p>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</p></td>
</tr> </tr>
<tr> <tr>
<th>year_min</th> <th>year_min</th>
@ -269,7 +277,7 @@
</tr> </tr>
<tr> <tr>
<th>model</th> <th>model</th>
<td><p>the statistical model of choice. Valid values are <code>"binomial"</code> (or <code>"binom"</code> or <code>"logit"</code>) or <code>"loglin"</code> or <code>"linear"</code> (or <code>"lin"</code>).</p></td> <td><p>the statistical model of choice. Valid values are <code>"binomial"</code> (or <code>"binom"</code> or <code>"logit"</code>) or <code>"loglin"</code> (or <code>"poisson"</code>) or <code>"linear"</code> (or <code>"lin"</code>).</p></td>
</tr> </tr>
<tr> <tr>
<th>I_as_R</th> <th>I_as_R</th>
@ -283,11 +291,21 @@
<th>info</th> <th>info</th>
<td><p>a logical to indicate whether textual analysis should be printed with the name and <code><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></code> of the statistical model.</p></td> <td><p>a logical to indicate whether textual analysis should be printed with the name and <code><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></code> of the statistical model.</p></td>
</tr> </tr>
<tr>
<th>x</th>
<td><p>the coordinates of points in the plot. Alternatively, a
single plotting structure, function or <em>any <span style="R">R</span> object with a
<code>plot</code> method</em> can be provided.</p></td>
</tr>
<tr>
<th>...</th>
<td><p>parameters passed on to the <code>first_isolate</code> function</p></td>
</tr>
</table> </table>
<h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2> <h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2>
<p><code>data.frame</code> with columns:</p><ul> <p><code>data.frame</code> with extra class <code>"resistance_predict"</code> with columns:</p><ul>
<li><p><code>year</code></p></li> <li><p><code>year</code></p></li>
<li><p><code>value</code>, the same as <code>estimated</code> when <code>preserve_measurements = FALSE</code>, and a combination of <code>observed</code> and <code>estimated</code> otherwise</p></li> <li><p><code>value</code>, the same as <code>estimated</code> when <code>preserve_measurements = FALSE</code>, and a combination of <code>observed</code> and <code>estimated</code> otherwise</p></li>
<li><p><code>se_min</code>, the lower bound of the standard error with a minimum of <code>0</code> (so the standard error will never go below 0%)</p></li> <li><p><code>se_min</code>, the lower bound of the standard error with a minimum of <code>0</code> (so the standard error will never go below 0%)</p></li>
@ -311,37 +329,20 @@ On our website <a href='https://msberends.gitlab.io/AMR'>https://msberends.gitla
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span> <pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'># use it with base R:</span> <span class='no'>x</span> <span class='kw'>&lt;-</span> <span class='fu'>resistance_predict</span>(<span class='no'>septic_patients</span>, <span class='kw'>col_ab</span> <span class='kw'>=</span> <span class='st'>"amox"</span>, <span class='kw'>year_min</span> <span class='kw'>=</span> <span class='fl'>2010</span>)
<span class='fu'>resistance_predict</span>(<span class='kw'>tbl</span> <span class='kw'>=</span> <span class='no'>tbl</span>[<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/which'>which</a></span>(<span class='no'>first_isolate</span> <span class='kw'>==</span> <span class='fl'>TRUE</span> <span class='kw'>&amp;</span> <span class='no'>genus</span> <span class='kw'>==</span> <span class='st'>"Haemophilus"</span>),], <span class='fu'><a href='https://www.rdocumentation.org/packages/graphics/topics/plot'>plot</a></span>(<span class='no'>x</span>)
<span class='kw'>col_ab</span> <span class='kw'>=</span> <span class='st'>"amcl"</span>, <span class='kw'>col_date</span> <span class='kw'>=</span> <span class='st'>"date"</span>) <span class='fu'>ggplot_rsi_predict</span>(<span class='no'>x</span>)
<span class='co'># or use dplyr so you can actually read it:</span> <span class='co'># use dplyr so you can actually read it:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>dplyr</span>) <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>dplyr</span>)
<span class='no'>tbl</span> <span class='kw'>%&gt;%</span> <span class='no'>x</span> <span class='kw'>&lt;-</span> <span class='no'>septic_patients</span> <span class='kw'>%&gt;%</span>
<span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span>(<span class='no'>first_isolate</span> <span class='kw'>==</span> <span class='fl'>TRUE</span>, <span class='fu'><a href='first_isolate.html'>filter_first_isolate</a></span>() <span class='kw'>%&gt;%</span>
<span class='no'>genus</span> <span class='kw'>==</span> <span class='st'>"Haemophilus"</span>) <span class='kw'>%&gt;%</span> <span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span>(<span class='fu'><a href='mo_property.html'>mo_genus</a></span>(<span class='no'>mo</span>) <span class='kw'>==</span> <span class='st'>"Staphylococcus"</span>) <span class='kw'>%&gt;%</span>
<span class='fu'>resistance_predict</span>(<span class='no'>amcl</span>, <span class='no'>date</span>) <span class='fu'>resistance_predict</span>(<span class='st'>"peni"</span>)
<span class='co'># }</span><span class='co'># NOT RUN {</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/graphics/topics/plot'>plot</a></span>(<span class='no'>x</span>)
<span class='co'># real live example:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>dplyr</span>)
<span class='no'>septic_patients</span> <span class='kw'>%&gt;%</span>
<span class='co'># get bacteria properties like genus and species</span>
<span class='fu'><a href='join.html'>left_join_microorganisms</a></span>(<span class='st'>"mo"</span>) <span class='kw'>%&gt;%</span>
<span class='co'># calculate first isolates</span>
<span class='fu'><a href='https://dplyr.tidyverse.org/reference/mutate.html'>mutate</a></span>(<span class='kw'>first_isolate</span> <span class='kw'>=</span> <span class='fu'><a href='first_isolate.html'>first_isolate</a></span>(<span class='no'>.</span>)) <span class='kw'>%&gt;%</span>
<span class='co'># filter on first E. coli isolates</span>
<span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span>(<span class='no'>genus</span> <span class='kw'>==</span> <span class='st'>"Escherichia"</span>,
<span class='no'>species</span> <span class='kw'>==</span> <span class='st'>"coli"</span>,
<span class='no'>first_isolate</span> <span class='kw'>==</span> <span class='fl'>TRUE</span>) <span class='kw'>%&gt;%</span>
<span class='co'># predict resistance of cefotaxime for next years</span>
<span class='fu'>resistance_predict</span>(<span class='kw'>col_ab</span> <span class='kw'>=</span> <span class='st'>"cfot"</span>,
<span class='kw'>col_date</span> <span class='kw'>=</span> <span class='st'>"date"</span>,
<span class='kw'>year_max</span> <span class='kw'>=</span> <span class='fl'>2025</span>,
<span class='kw'>preserve_measurements</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>,
<span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>0</span>)
<span class='co'># create nice plots with ggplot</span> <span class='co'># create nice plots with ggplot yourself</span>
<span class='kw'>if</span> (!<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>require</a></span>(<span class='no'>ggplot2</span>)) { <span class='kw'>if</span> (!<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>require</a></span>(<span class='no'>ggplot2</span>)) {
<span class='no'>data</span> <span class='kw'>&lt;-</span> <span class='no'>septic_patients</span> <span class='kw'>%&gt;%</span> <span class='no'>data</span> <span class='kw'>&lt;-</span> <span class='no'>septic_patients</span> <span class='kw'>%&gt;%</span>
@ -349,7 +350,7 @@ On our website <a href='https://msberends.gitlab.io/AMR'>https://msberends.gitla
<span class='fu'>resistance_predict</span>(<span class='kw'>col_ab</span> <span class='kw'>=</span> <span class='st'>"amox"</span>, <span class='fu'>resistance_predict</span>(<span class='kw'>col_ab</span> <span class='kw'>=</span> <span class='st'>"amox"</span>,
<span class='kw'>col_date</span> <span class='kw'>=</span> <span class='st'>"date"</span>, <span class='kw'>col_date</span> <span class='kw'>=</span> <span class='st'>"date"</span>,
<span class='kw'>info</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='kw'>info</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>15</span>) <span class='kw'>minimum</span> <span class='kw'>=</span> <span class='fl'>15</span>)
<span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span>(<span class='no'>data</span>, <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span>(<span class='no'>data</span>,
<span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span>(<span class='kw'>x</span> <span class='kw'>=</span> <span class='no'>year</span>)) + <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span>(<span class='kw'>x</span> <span class='kw'>=</span> <span class='no'>year</span>)) +

View File

@ -26,7 +26,7 @@ filter_first_weighted_isolate(tbl, col_date = NULL,
\arguments{ \arguments{
\item{tbl}{a \code{data.frame} containing isolates.} \item{tbl}{a \code{data.frame} containing isolates.}
\item{col_date}{column name of the result date (or date that is was received on the lab), defaults to the first column of class \code{Date}} \item{col_date}{column name of the result date (or date that is was received on the lab), defaults to the first column of with a date class}
\item{col_patient_id}{column name of the unique IDs of the patients, defaults to the first column that starts with 'patient' or 'patid' (case insensitive)} \item{col_patient_id}{column name of the unique IDs of the patients, defaults to the first column that starts with 'patient' or 'patid' (case insensitive)}

View File

@ -3,23 +3,32 @@
\name{resistance_predict} \name{resistance_predict}
\alias{resistance_predict} \alias{resistance_predict}
\alias{rsi_predict} \alias{rsi_predict}
\alias{plot.resistance_predict}
\alias{ggplot_rsi_predict}
\title{Predict antimicrobial resistance} \title{Predict antimicrobial resistance}
\usage{ \usage{
resistance_predict(tbl, col_ab, col_date, year_min = NULL, resistance_predict(tbl, col_ab, col_date = NULL, year_min = NULL,
year_max = NULL, year_every = 1, minimum = 30, year_max = NULL, year_every = 1, minimum = 30,
model = "binomial", I_as_R = TRUE, preserve_measurements = TRUE, model = "binomial", I_as_R = TRUE, preserve_measurements = TRUE,
info = TRUE) info = TRUE)
rsi_predict(tbl, col_ab, col_date, year_min = NULL, year_max = NULL, rsi_predict(tbl, col_ab, col_date = NULL, year_min = NULL,
year_every = 1, minimum = 30, model = "binomial", I_as_R = TRUE, year_max = NULL, year_every = 1, minimum = 30,
preserve_measurements = TRUE, info = TRUE) model = "binomial", I_as_R = TRUE, preserve_measurements = TRUE,
info = TRUE)
\method{plot}{resistance_predict}(x,
main = paste("Resistance prediction of", attributes(x)$ab), ...)
ggplot_rsi_predict(x, main = paste("Resistance prediction of",
attributes(x)$ab), ...)
} }
\arguments{ \arguments{
\item{tbl}{a \code{data.frame} containing isolates.} \item{tbl}{a \code{data.frame} containing isolates.}
\item{col_ab}{column name of \code{tbl} with antimicrobial interpretations (\code{R}, \code{I} and \code{S})} \item{col_ab}{column name of \code{tbl} with antimicrobial interpretations (\code{R}, \code{I} and \code{S})}
\item{col_date}{column name of the date, will be used to calculate years if this column doesn't consist of years already} \item{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}
\item{year_min}{lowest year to use in the prediction model, dafaults to the lowest year in \code{col_date}} \item{year_min}{lowest year to use in the prediction model, dafaults to the lowest year in \code{col_date}}
@ -29,16 +38,22 @@ rsi_predict(tbl, col_ab, col_date, year_min = NULL, year_max = NULL,
\item{minimum}{minimal amount of available isolates per year to include. Years containing less observations will be estimated by the model.} \item{minimum}{minimal amount of available isolates per year to include. Years containing less observations will be estimated by the model.}
\item{model}{the statistical model of choice. Valid values are \code{"binomial"} (or \code{"binom"} or \code{"logit"}) or \code{"loglin"} or \code{"linear"} (or \code{"lin"}).} \item{model}{the statistical model of choice. Valid values are \code{"binomial"} (or \code{"binom"} or \code{"logit"}) or \code{"loglin"} (or \code{"poisson"}) or \code{"linear"} (or \code{"lin"}).}
\item{I_as_R}{a logical to indicate whether values \code{I} should be treated as \code{R}} \item{I_as_R}{a logical to indicate whether values \code{I} should be treated as \code{R}}
\item{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 \code{NA}.} \item{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 \code{NA}.}
\item{info}{a logical to indicate whether textual analysis should be printed with the name and \code{\link{summary}} of the statistical model.} \item{info}{a logical to indicate whether textual analysis should be printed with the name and \code{\link{summary}} of the statistical model.}
\item{x}{the coordinates of points in the plot. Alternatively, a
single plotting structure, function or \emph{any \R object with a
\code{plot} method} can be provided.}
\item{...}{parameters passed on to the \code{first_isolate} function}
} }
\value{ \value{
\code{data.frame} with columns: \code{data.frame} with extra class \code{"resistance_predict"} with columns:
\itemize{ \itemize{
\item{\code{year}} \item{\code{year}}
\item{\code{value}, the same as \code{estimated} when \code{preserve_measurements = FALSE}, and a combination of \code{observed} and \code{estimated} otherwise} \item{\code{value}, the same as \code{estimated} when \code{preserve_measurements = FALSE}, and a combination of \code{observed} and \code{estimated} otherwise}
@ -59,39 +74,20 @@ On our website \url{https://msberends.gitlab.io/AMR} you can find \href{https://
} }
\examples{ \examples{
\dontrun{ x <- resistance_predict(septic_patients, col_ab = "amox", year_min = 2010)
# use it with base R: plot(x)
resistance_predict(tbl = tbl[which(first_isolate == TRUE & genus == "Haemophilus"),], ggplot_rsi_predict(x)
col_ab = "amcl", col_date = "date")
# or use dplyr so you can actually read it: # use dplyr so you can actually read it:
library(dplyr) library(dplyr)
tbl \%>\% x <- septic_patients \%>\%
filter(first_isolate == TRUE, filter_first_isolate() \%>\%
genus == "Haemophilus") \%>\% filter(mo_genus(mo) == "Staphylococcus") \%>\%
resistance_predict(amcl, date) resistance_predict("peni")
} plot(x)
# real live example: # create nice plots with ggplot yourself
library(dplyr)
septic_patients \%>\%
# get bacteria properties like genus and species
left_join_microorganisms("mo") \%>\%
# calculate first isolates
mutate(first_isolate = first_isolate(.)) \%>\%
# filter on first E. coli isolates
filter(genus == "Escherichia",
species == "coli",
first_isolate == TRUE) \%>\%
# predict resistance of cefotaxime for next years
resistance_predict(col_ab = "cfot",
col_date = "date",
year_max = 2025,
preserve_measurements = TRUE,
minimum = 0)
# create nice plots with ggplot
if (!require(ggplot2)) { if (!require(ggplot2)) {
data <- septic_patients \%>\% data <- septic_patients \%>\%
@ -99,7 +95,7 @@ if (!require(ggplot2)) {
resistance_predict(col_ab = "amox", resistance_predict(col_ab = "amox",
col_date = "date", col_date = "date",
info = FALSE, info = FALSE,
minimum = 15) minimum = 15)
ggplot(data, ggplot(data,
aes(x = year)) + aes(x = year)) +

View File

@ -32,6 +32,11 @@ test_that("prediction of rsi works", {
# amox resistance will increase according to data set `septic_patients` # amox resistance will increase according to data set `septic_patients`
expect_true(amox_R[3] < amox_R[20]) expect_true(amox_R[3] < amox_R[20])
x <- resistance_predict(septic_patients, col_ab = "amox", year_min = 2010)
plot(x)
ggplot_rsi_predict(x)
expect_error(ggplot_rsi_predict(septic_patients))
library(dplyr) library(dplyr)
expect_output(rsi_predict(tbl = filter(septic_patients, mo == "B_ESCHR_COL"), expect_output(rsi_predict(tbl = filter(septic_patients, mo == "B_ESCHR_COL"),