1
0
mirror of https://github.com/msberends/AMR.git synced 2025-07-09 06:02:01 +02:00

update dependencies

This commit is contained in:
2018-04-02 16:05:09 +02:00
parent cee64ef050
commit 07bdd61241
15 changed files with 172 additions and 221 deletions

View File

@ -48,7 +48,7 @@
#' cfur = "-", # Cefuroxime
#' stringsAsFactors = FALSE)
#' a
#'
#'
#' b <- EUCAST_rules(a)
#' b
EUCAST_rules <- function(tbl,
@ -114,71 +114,11 @@ EUCAST_rules <- function(tbl,
vanc = 'vanc') {
EUCAST_VERSION <- "3.1"
# support using columns as objects; the tidyverse way
amcl <- quasiquotate(deparse(substitute(amcl)), amcl)
amik <- quasiquotate(deparse(substitute(amik)), amik)
amox <- quasiquotate(deparse(substitute(amox)), amox)
ampi <- quasiquotate(deparse(substitute(ampi)), ampi)
azit <- quasiquotate(deparse(substitute(azit)), azit)
aztr <- quasiquotate(deparse(substitute(aztr)), aztr)
cefa <- quasiquotate(deparse(substitute(cefa)), cefa)
cfra <- quasiquotate(deparse(substitute(cfra)), cfra)
cfep <- quasiquotate(deparse(substitute(cfep)), cfep)
cfot <- quasiquotate(deparse(substitute(cfot)), cfot)
cfox <- quasiquotate(deparse(substitute(cfox)), cfox)
cfta <- quasiquotate(deparse(substitute(cfta)), cfta)
cftr <- quasiquotate(deparse(substitute(cftr)), cftr)
cfur <- quasiquotate(deparse(substitute(cfur)), cfur)
chlo <- quasiquotate(deparse(substitute(chlo)), chlo)
cipr <- quasiquotate(deparse(substitute(cipr)), cipr)
clar <- quasiquotate(deparse(substitute(clar)), clar)
clin <- quasiquotate(deparse(substitute(clin)), clin)
clox <- quasiquotate(deparse(substitute(clox)), clox)
coli <- quasiquotate(deparse(substitute(coli)), coli)
czol <- quasiquotate(deparse(substitute(czol)), czol)
dapt <- quasiquotate(deparse(substitute(dapt)), dapt)
doxy <- quasiquotate(deparse(substitute(doxy)), doxy)
erta <- quasiquotate(deparse(substitute(erta)), erta)
eryt <- quasiquotate(deparse(substitute(eryt)), eryt)
fosf <- quasiquotate(deparse(substitute(fosf)), fosf)
fusi <- quasiquotate(deparse(substitute(fusi)), fusi)
gent <- quasiquotate(deparse(substitute(gent)), gent)
imip <- quasiquotate(deparse(substitute(imip)), imip)
kana <- quasiquotate(deparse(substitute(kana)), kana)
levo <- quasiquotate(deparse(substitute(levo)), levo)
linc <- quasiquotate(deparse(substitute(linc)), linc)
line <- quasiquotate(deparse(substitute(line)), line)
mero <- quasiquotate(deparse(substitute(mero)), mero)
mino <- quasiquotate(deparse(substitute(mino)), mino)
moxi <- quasiquotate(deparse(substitute(moxi)), moxi)
nali <- quasiquotate(deparse(substitute(nali)), nali)
neom <- quasiquotate(deparse(substitute(neom)), neom)
neti <- quasiquotate(deparse(substitute(neti)), neti)
nitr <- quasiquotate(deparse(substitute(nitr)), nitr)
novo <- quasiquotate(deparse(substitute(novo)), novo)
norf <- quasiquotate(deparse(substitute(norf)), norf)
oflo <- quasiquotate(deparse(substitute(oflo)), oflo)
peni <- quasiquotate(deparse(substitute(peni)), peni)
pita <- quasiquotate(deparse(substitute(pita)), pita)
poly <- quasiquotate(deparse(substitute(poly)), poly)
qida <- quasiquotate(deparse(substitute(qida)), qida)
rifa <- quasiquotate(deparse(substitute(rifa)), rifa)
roxi <- quasiquotate(deparse(substitute(roxi)), roxi)
siso <- quasiquotate(deparse(substitute(siso)), siso)
teic <- quasiquotate(deparse(substitute(teic)), teic)
tetr <- quasiquotate(deparse(substitute(tetr)), tetr)
tica <- quasiquotate(deparse(substitute(tica)), tica)
tige <- quasiquotate(deparse(substitute(tige)), tige)
tobr <- quasiquotate(deparse(substitute(tobr)), tobr)
trim <- quasiquotate(deparse(substitute(trim)), trim)
trsu <- quasiquotate(deparse(substitute(trsu)), trsu)
vanc <- quasiquotate(deparse(substitute(vanc)), vanc)
if (!col_bactid %in% colnames(tbl)) {
stop('Column ', col_bactid, ' not found.')
}
# check columns
col.list <- c(amcl, amik, amox, ampi, azit, aztr, cefa, cfra, cfep, cfot,
cfox, cfta, cftr, cfur, chlo, cipr, clar, clin, clox, coli,
@ -209,7 +149,7 @@ EUCAST_rules <- function(tbl,
call. = FALSE)
}
}
amcl <- col.list[1]
amik <- col.list[2]
amox <- col.list[3]
@ -268,10 +208,10 @@ EUCAST_rules <- function(tbl,
trim <- col.list[56]
trsu <- col.list[57]
vanc <- col.list[58]
total <- 0
total_rows <- integer(0)
# helper function for editing the table
edit_rsi <- function(to, rows, cols) {
cols <- cols[!is.na(cols)]
@ -281,12 +221,12 @@ EUCAST_rules <- function(tbl,
total_rows <<- c(total_rows, rows)
}
}
# join to microorganisms table
joinby <- colnames(AMR::microorganisms)[1]
names(joinby) <- col_bactid
tbl <- tbl %>% left_join(y = AMR::microorganisms, by = joinby, suffix = c("_tempmicroorganisms", ""))
# antibiotic classes
aminoglycosides <- c(tobr, gent, kana, neom, neti, siso)
tetracyclines <- c(doxy, mino, tetr) # since EUCAST v3.1 tige(cycline) is set apart
@ -299,7 +239,7 @@ EUCAST_rules <- function(tbl,
aminopenicillins <- c(ampi, amox)
ureidopenicillins <- pita # should officially also be azlo and mezlo
fluoroquinolones <- c(oflo, cipr, norf, levo, moxi)
if (info == TRUE) {
cat(
paste0(
@ -308,7 +248,7 @@ EUCAST_rules <- function(tbl,
' rows according to "EUCAST Expert Rules Version ', EUCAST_VERSION, '"\n')
)
}
# Table 1: Intrinsic resistance in Enterobacteriaceae ----
if (info == TRUE) {
cat('...Table 1: Intrinsic resistance in Enterobacteriaceae\n')
@ -378,8 +318,8 @@ EUCAST_rules <- function(tbl,
edit_rsi(to = 'R',
rows = which(tbl$fullname %like% '^Yersinia pseudotuberculosis'),
cols = c(poly, coli))
# Table 2: Intrinsic resistance in non-fermentative Gram-negative bacteria ----
if (info == TRUE) {
cat('...Table 2: Intrinsic resistance in non-fermentative Gram-negative bacteria\n')
@ -426,8 +366,8 @@ EUCAST_rules <- function(tbl,
edit_rsi(to = 'R',
rows = which(tbl$fullname %like% '^Stenotrophomonas maltophilia'),
cols = c(aminopenicillins, amcl, tica, pita, czol, cfot, cftr, cfta, aztr, erta, imip, mero, aminoglycosides, trim, fosf, tetr))
# Table 3: Intrinsic resistance in other Gram-negative bacteria ----
if (info == TRUE) {
cat('...Table 3: Intrinsic resistance in other Gram-negative bacteria\n')
@ -458,8 +398,8 @@ EUCAST_rules <- function(tbl,
edit_rsi(to = 'R',
rows = which(tbl$fullname %like% '^Campylobacter (jejuni|coli)'),
cols = c(fusi, streptogramins, trim))
# Table 4: Intrinsic resistance in Gram-positive bacteria ----
if (info == TRUE) {
cat('...Table 4: Intrinsic resistance in Gram-positive bacteria\n')
@ -513,7 +453,7 @@ EUCAST_rules <- function(tbl,
edit_rsi(to = 'R',
rows = which(tbl$fullname %like% '^Clostridium (ramosum|innocuum)'),
cols = vanc)
# Table 8: Interpretive rules for B-lactam agents and Gram-positive cocci ----
if (info == TRUE) {
cat('...Table 8: Interpretive rules for B-lactam agents and Gram-positive cocci\n')
@ -538,7 +478,7 @@ EUCAST_rules <- function(tbl,
& tbl[, amox] == 'R'),
cols = c(ureidopenicillins, carbapenems))
}
# Table 9: Interpretive rules for B-lactam agents and Gram-negative rods ----
if (info == TRUE) {
cat('...Table 9: Interpretive rules for B-lactam agents and Gram-negative rods\n')
@ -551,7 +491,7 @@ EUCAST_rules <- function(tbl,
& tbl[, pita] == 'S'),
cols = pita)
}
# Table 10: Interpretive rules for B-lactam agents and other Gram-negative bacteria ----
if (info == TRUE) {
cat('...Table 10: Interpretive rules for B-lactam agents and other Gram-negative bacteria\n')
@ -564,7 +504,7 @@ EUCAST_rules <- function(tbl,
# & tbl[, ampi] == 'R'),
# cols = c(ampi, amox, amcl, pita, cfur))
}
# Table 11: Interpretive rules for macrolides, lincosamides, and streptogramins ----
if (info == TRUE) {
cat('...Table 11: Interpretive rules for macrolides, lincosamides, and streptogramins\n')
@ -578,7 +518,7 @@ EUCAST_rules <- function(tbl,
tbl[, clar] <- tbl[, eryt]
}
}
# Table 12: Interpretive rules for aminoglycosides ----
if (info == TRUE) {
cat('...Table 12: Interpretive rules for aminoglycosides\n')
@ -613,8 +553,8 @@ EUCAST_rules <- function(tbl,
& tbl[, gent] == 'R'),
cols = tobr)
}
# Table 13: Interpretive rules for quinolones ----
if (info == TRUE) {
cat('...Table 13: Interpretive rules for quinolones\n')
@ -647,8 +587,8 @@ EUCAST_rules <- function(tbl,
& tbl[, cipr] == 'R'),
cols = fluoroquinolones)
}
# Other ----
if (info == TRUE) {
cat('...Non-EUCAST: trim = R where trsu = R and ampi = R where amcl = R\n')
@ -666,21 +606,21 @@ EUCAST_rules <- function(tbl,
if (!is.na(ampi) & !is.na(amox)) {
tbl[, amox] <- tbl %>% pull(ampi)
}
# Remove added columns again
microorganisms.ncol <- ncol(AMR::microorganisms) - 2
tbl.ncol <- ncol(tbl)
tbl <- tbl %>% select(-c((tbl.ncol - microorganisms.ncol):tbl.ncol))
# and remove added suffices
colnames(tbl) <- gsub("_tempmicroorganisms", "", colnames(tbl))
if (info == TRUE) {
cat('Done.\n\nEUCAST Expert rules applied to',
total_rows %>% unique() %>% length() %>% format(big.mark = ","),
'different rows (isolates); edited a total of',
total %>% format(big.mark = ","), 'test results.\n\n')
}
tbl
}
@ -698,12 +638,12 @@ interpretive_reading <- function(...) {
#' @importFrom dplyr %>% filter select
#' @seealso \code{\link{microorganisms}}
mo_property <- function(bactid, property = 'fullname') {
mocode <- as.character(bactid)
for (i in 1:length(mocode)) {
bug <- mocode[i]
if (!is.na(bug)) {
result = tryCatch({
mocode[i] <-
@ -720,7 +660,7 @@ mo_property <- function(bactid, property = 'fullname') {
}
})
}
}
mocode
}

View File

@ -29,18 +29,23 @@
#' rsi_data <- as.rsi(c(rep("S", 474), rep("I", 36), rep("R", 370)))
#' rsi_data <- as.rsi(c(rep("S", 474), rep("I", 36), rep("R", 370), "A", "B", "C"))
#' is.rsi(rsi_data)
#'
#'
#' plot(rsi_data) # for percentages
#' barplot(rsi_data) # for frequencies
as.rsi <- function(x) {
if (is.rsi(x)) {
x
} else {
x <- x %>% unlist()
x.bak <- x
na_before <- x[is.na(x) | x == ''] %>% length()
# remove all spaces
x <- gsub(' {2,55}', '', x)
# disallow more than 3 characters
x[nchar(x) > 3] <- NA
# remove all invalid characters
x <- gsub('[^RSI]+', '', x %>% toupper())
# needed for UMCG in cases of "S;S" but also "S;I"; the latter will be NA:
x <- gsub('^S+$', 'S', x)
@ -48,7 +53,7 @@ as.rsi <- function(x) {
x <- gsub('^R+$', 'R', x)
x[!x %in% c('S', 'I', 'R')] <- NA
na_after <- x[is.na(x) | x == ''] %>% length()
if (na_before != na_after) {
list_missing <- x.bak[is.na(x) & !is.na(x.bak) & x.bak != ''] %>%
unique() %>%
@ -59,7 +64,7 @@ as.rsi <- function(x) {
'%) that were invalid antimicrobial interpretations: ',
list_missing, call. = FALSE)
}
x <- x %>% toupper() %>% factor(levels = c("S", "I", "R"), ordered = TRUE)
class(x) <- c('rsi', 'ordered', 'factor')
attr(x, 'package') <- 'AMR'
@ -128,7 +133,7 @@ summary.rsi <- function(object, ...) {
#' @noRd
plot.rsi <- function(x, ...) {
x_name <- deparse(substitute(x))
data <- data.frame(x = x,
y = 1,
stringsAsFactors = TRUE) %>%
@ -137,7 +142,7 @@ plot.rsi <- function(x, ...) {
filter(!is.na(x)) %>%
mutate(s = round((n / sum(n)) * 100, 1))
data$x <- factor(data$x, levels = c('S', 'I', 'R'), ordered = TRUE)
ymax <- if_else(max(data$s) > 95, 105, 100)
plot(x = data$x,
@ -154,7 +159,7 @@ plot.rsi <- function(x, ...) {
axis(side = 1, at = 1:n_distinct(data$x), labels = levels(data$x), lwd = 0)
# y axis, 0-100%
axis(side = 2, at = seq(0, 100, 5))
text(x = data$x,
y = data$s + 4,
labels = paste0(data$s, '% (n = ', data$n, ')'))
@ -169,7 +174,7 @@ plot.rsi <- function(x, ...) {
barplot.rsi <- function(height, ...) {
x <- height
x_name <- deparse(substitute(height))
data <- data.frame(rsi = x, cnt = 1) %>%
group_by(rsi) %>%
summarise(cnt = sum(cnt)) %>%
@ -199,7 +204,7 @@ barplot.rsi <- function(height, ...) {
#' @examples
#' mic_data <- as.mic(c(">=32", "1.0", "1", "1.00", 8, "<=0.128", "8", "16", "16"))
#' is.mic(mic_data)
#'
#'
#' plot(mic_data)
#' barplot(mic_data)
as.mic <- function(x, na.rm = FALSE) {
@ -211,7 +216,7 @@ as.mic <- function(x, na.rm = FALSE) {
x <- x[!is.na(x)]
}
x.bak <- x
# comma to dot
x <- gsub(',', '.', x, fixed = TRUE)
# starting dots must start with 0
@ -224,7 +229,7 @@ as.mic <- function(x, na.rm = FALSE) {
x <- gsub('[^0-9]$', '', x)
# remove last zeroes
x <- gsub('[.]?0+$', '', x)
lvls <- c("<0.002", "<=0.002", "0.002", ">=0.002", ">0.002",
"<0.003", "<=0.003", "0.003", ">=0.003", ">0.003",
"<0.004", "<=0.004", "0.004", ">=0.004", ">0.004",
@ -282,11 +287,11 @@ as.mic <- function(x, na.rm = FALSE) {
"<512", "<=512", "512", ">=512", ">512",
"<1024", "<=1024", "1024", ">=1024", ">1024")
x <- x %>% as.character()
na_before <- x[is.na(x) | x == ''] %>% length()
x[!x %in% lvls] <- NA
na_after <- x[is.na(x) | x == ''] %>% length()
if (na_before != na_after) {
list_missing <- x.bak[is.na(x) & !is.na(x.bak) & x.bak != ''] %>%
unique() %>%
@ -297,7 +302,7 @@ as.mic <- function(x, na.rm = FALSE) {
'%) that were invalid MICs: ',
list_missing, call. = FALSE)
}
x <- factor(x = x,
levels = lvls,
ordered = TRUE)
@ -407,7 +412,7 @@ create_barplot_mic <- function(x, x_name, ...) {
barplot(table(droplevels(x)),
ylab = 'Frequency',
xlab = 'MIC value',
main = paste('MIC values of', x_name),
main = paste('MIC values of', x_name),
axes = FALSE,
...)
axis(2, seq(0, max(data$cnt)))

View File

@ -31,22 +31,3 @@
percent <- function(x, round = 1, ...) {
base::paste0(base::round(x * 100, digits = round), "%")
}
# No export, no Rd
quasiquotate <- function(deparsed, parsed) {
# when text: remove first and last "
if (any(deparsed %like% '^".+"$' | deparsed %like% "^'.+'$")) {
deparsed <- deparsed %>% substr(2, nchar(.) - 1)
}
# apply if needed
if (any(!deparsed %like% '[[$:()]'
& !deparsed %in% c('""', "''", "", # empty text
".", ".data", # dplyr references
"TRUE", "FALSE", # logicals
"NA", "NaN", "NULL", # empty values
ls(.GlobalEnv)))) {
deparsed
} else {
parsed
}
}

View File

@ -41,7 +41,7 @@
#' library(dplyr)
#' # calculate current empiric therapy of Helicobacter gastritis:
#' my_table %>%
#' filter(first_isolate == TRUE,
#' filter(first_isolate == TRUE,
#' genus == "Helicobacter") %>%
#' rsi_df(ab = c("amox", "metr"))
#' }
@ -55,7 +55,7 @@ rsi_df <- function(tbl,
# in case tbl$interpretation already exists:
interpretations_to_check <- paste(interpretation, collapse = "")
# validate:
if (min(grepl('^[a-z]{3,4}$', ab)) == 0 &
min(grepl('^rsi[1-2]$', ab)) == 0) {
@ -71,7 +71,7 @@ rsi_df <- function(tbl,
warning('Dataset contains isolates from the Intensive Care. Exclude them from proper epidemiological analysis.')
}
}
# transform when checking for different results
if (interpretations_to_check %in% c('SI', 'IS')) {
for (i in 1:length(ab)) {
@ -101,7 +101,7 @@ rsi_df <- function(tbl,
denominator <- tbl %>%
filter(pull(., ab[1]) %in% c("S", "I", "R")) %>%
nrow()
} else if (length(ab) == 2) {
numerator <- tbl %>%
filter_at(vars(ab[1], ab[2]),
@ -109,12 +109,12 @@ rsi_df <- function(tbl,
filter_at(vars(ab[1], ab[2]),
all_vars(. %in% c("S", "R", "I"))) %>%
nrow()
denominator <- tbl %>%
filter_at(vars(ab[1], ab[2]),
all_vars(. %in% c("S", "R", "I"))) %>%
nrow()
} else if (length(ab) == 3) {
numerator <- tbl %>%
filter_at(vars(ab[1], ab[2], ab[3]),
@ -122,16 +122,16 @@ rsi_df <- function(tbl,
filter_at(vars(ab[1], ab[2], ab[3]),
all_vars(. %in% c("S", "R", "I"))) %>%
nrow()
denominator <- tbl %>%
filter_at(vars(ab[1], ab[2], ab[3]),
all_vars(. %in% c("S", "R", "I"))) %>%
nrow()
} else {
stop('Maximum of 3 drugs allowed.')
}
# build text part
if (info == TRUE) {
cat('n =', denominator)
@ -147,7 +147,7 @@ rsi_df <- function(tbl,
info.txt2 <- gsub('rsi1', 'this drug', info.txt2, fixed = TRUE)
cat(paste0(' (of ', nrow(tbl), ' in total; ', info.txt1, ' tested on ', info.txt2, ')\n'))
}
# calculate and format
y <- numerator / denominator
if (percent == TRUE) {
@ -159,7 +159,7 @@ rsi_df <- function(tbl,
}
y <- NA
}
# output
y
}
@ -178,14 +178,14 @@ rsi_df <- function(tbl,
#' tbl %>%
#' group_by(hospital) %>%
#' summarise(cipr = rsi(cipr))
#'
#'
#' tbl %>%
#' group_by(year, hospital) %>%
#' summarise(
#' isolates = n(),
#' cipro = rsi(cipr %>% as.rsi(), percent = TRUE),
#' amoxi = rsi(amox %>% as.rsi(), percent = TRUE))
#'
#'
#' rsi(as.rsi(isolates$amox))
#'
#' rsi(as.rsi(isolates$amcl), interpretation = "S")
@ -207,12 +207,12 @@ rsi <- function(ab1, ab2 = NA, interpretation = 'IR', minimum = 30, percent = FA
if (!ab2.name %like% '^[a-z]{3,4}$') {
ab2.name <- 'rsi2'
}
interpretation <- paste(interpretation, collapse = "")
tbl <- tibble(rsi1 = ab1, rsi2 = ab2)
colnames(tbl) <- c(ab1.name, ab2.name)
if (length(ab2) == 1) {
return(rsi_df(tbl = tbl,
ab = ab1.name,
@ -260,7 +260,7 @@ rsi <- function(ab1, ab2 = NA, interpretation = 'IR', minimum = 30, percent = FA
#' # use it directly:
#' rsi_predict(tbl = tbl[which(first_isolate == TRUE & genus == "Haemophilus"),],
#' col_ab = "amcl", col_date = "date")
#'
#'
#' # or with dplyr so you can actually read it:
#' library(dplyr)
#' tbl %>%
@ -274,22 +274,22 @@ rsi <- function(ab1, ab2 = NA, interpretation = 'IR', minimum = 30, percent = FA
#' library(dplyr)
#' septic_patients %>%
#' # get bacteria properties like genus and species
#' left_join_microorganisms("bactid") %>%
#' left_join_microorganisms("bactid") %>%
#' # calculate first isolates
#' mutate(first_isolate =
#' mutate(first_isolate =
#' first_isolate(.,
#' "date",
#' "patient_id",
#' "bactid",
#' col_specimen = NA,
#' col_icu = NA)) %>%
#' col_icu = NA)) %>%
#' # filter on first E. coli isolates
#' filter(genus == "Escherichia",
#' species == "coli",
#' filter(genus == "Escherichia",
#' species == "coli",
#' first_isolate == TRUE) %>%
#' # predict resistance of cefotaxime for next years
#' rsi_predict(col_ab = cfot,
#' col_date = date,
#' rsi_predict(col_ab = "cfot",
#' col_date = "date",
#' year_max = 2025,
#' preserve_measurements = FALSE)
#'
@ -302,16 +302,15 @@ rsi_predict <- function(tbl,
I_as_R = TRUE,
preserve_measurements = TRUE,
info = TRUE) {
if (nrow(tbl) == 0) {
stop('This table does not contain any observations.')
}
col_ab <- quasiquotate(deparse(substitute(col_ab)), col_ab)
if (!col_ab %in% colnames(tbl)) {
stop('Column ', col_ab, ' not found.')
}
col_date <- quasiquotate(deparse(substitute(col_date)), col_date)
if (!col_date %in% colnames(tbl)) {
stop('Column ', col_date, ' not found.')
}
@ -327,7 +326,7 @@ rsi_predict <- function(tbl,
if (!all(tbl %>% pull(col_ab) %>% as.rsi() %in% c(NA, 'S', 'I', 'R'))) {
stop('Column ', col_ab, ' must contain antimicrobial interpretations (S, I, R).')
}
year <- function(x) {
if (all(grepl('^[0-9]{4}$', x))) {
x
@ -335,9 +334,9 @@ rsi_predict <- function(tbl,
as.integer(format(as.Date(x), '%Y'))
}
}
years_predict <- seq(from = min(year(tbl %>% pull(col_date))), to = year_max, by = year_every)
df <- tbl %>%
mutate(year = year(tbl %>% pull(col_date))) %>%
group_by_at(c('year', col_ab)) %>%
@ -345,7 +344,7 @@ rsi_predict <- function(tbl,
colnames(df) <- c('year', 'antibiotic', 'count')
df <- df %>%
reshape2::dcast(year ~ antibiotic, value.var = 'count')
if (model %in% c('binomial', 'binom', 'logit')) {
logitmodel <- with(df, glm(cbind(R, S) ~ year, family = binomial))
if (info == TRUE) {
@ -353,11 +352,11 @@ rsi_predict <- function(tbl,
cat('\n------------------------------------------------------------\n')
print(summary(logitmodel))
}
predictmodel <- stats::predict(logitmodel, newdata = with(df, list(year = years_predict)), type = "response", se.fit = TRUE)
prediction <- predictmodel$fit
se <- predictmodel$se.fit
} else if (model == 'loglin') {
loglinmodel <- with(df, glm(R ~ year, family = poisson))
if (info == TRUE) {
@ -365,11 +364,11 @@ rsi_predict <- function(tbl,
cat('\n--------------------------------------------------------------\n')
print(summary(loglinmodel))
}
predictmodel <- stats::predict(loglinmodel, newdata = with(df, list(year = years_predict)), type = "response", se.fit = TRUE)
prediction <- predictmodel$fit
se <- predictmodel$se.fit
} else if (model %in% c('lin', 'linear')) {
linmodel <- with(df, lm((R / (R + S)) ~ year))
if (info == TRUE) {
@ -377,36 +376,36 @@ rsi_predict <- function(tbl,
cat('\n-----------------------\n')
print(summary(linmodel))
}
predictmodel <- stats::predict(linmodel, newdata = with(df, list(year = years_predict)), se.fit = TRUE)
prediction <- predictmodel$fit
se <- predictmodel$se.fit
} else {
stop('No valid model selected.')
}
# prepare the output dataframe
prediction <- data.frame(year = years_predict, probR = prediction, stringsAsFactors = FALSE)
prediction$se_min <- prediction$probR - se
prediction$se_max <- prediction$probR + se
if (model == 'loglin') {
prediction$probR <- prediction$probR %>%
format(scientific = FALSE) %>%
as.integer()
prediction$se_min <- prediction$se_min %>% as.integer()
prediction$se_max <- prediction$se_max %>% as.integer()
colnames(prediction) <- c('year', 'amountR', 'se_max', 'se_min')
} else {
prediction$se_max[which(prediction$se_max > 1)] <- 1
}
prediction$se_min[which(prediction$se_min < 0)] <- 0
total <- prediction
if (preserve_measurements == TRUE) {
# geschatte data vervangen door gemeten data
if (I_as_R == TRUE) {
@ -424,10 +423,10 @@ rsi_predict <- function(tbl,
stringsAsFactors = FALSE)
colnames(measurements) <- colnames(prediction)
prediction <- prediction %>% filter(!year %in% df$year)
total <- rbind(measurements, prediction)
}
total
}