Introduction
@@ -217,21 +217,21 @@As with many uses in R, we need some additional packages for AMR analysis. Our package works closely together with the tidyverse packages dplyr and ggplot2 by Dr Hadley Wickham. The tidyverse tremendously improves the way we conduct data science - it allows for a very natural way of writing syntaxes and creating beautiful plots in R.
Our AMR package depends on these packages and even extends their use and functions.
library(dplyr)
-library(ggplot2)
-library(AMR)
-
-# (if not yet installed, install with:)
-# install.packages(c("tidyverse", "AMR"))library(dplyr)
+library(ggplot2)
+library(AMR)
+
+# (if not yet installed, install with:)
+# install.packages(c("tidyverse", "AMR"))Patients
To start with patients, we need a unique list of patients.
- +The LETTERS object is available in R - it’s a vector with 26 characters: A to Z. The patients object we just created is now a vector of length 260, with values (patient IDs) varying from A1 to Z10. Now we we also set the gender of our patients, by putting the ID and the gender in a table:
The first 135 patient IDs are now male, the other 125 are female.
Dates
Let’s pretend that our data consists of blood cultures isolates from between 1 January 2010 and 1 January 2018.
- +This dates object now contains all days in our date range.
Microorganisms
For this tutorial, we will uses four different microorganisms: Escherichia coli, Staphylococcus aureus, Streptococcus pneumoniae, and Klebsiella pneumoniae:
- +Other variables
For completeness, we can also add the hospital where the patients was admitted and we need to define valid antibmicrobial results for our randomisation:
- +Put everything together
Using the sample() function, we can randomly select items from all objects we defined earlier. To let our fake data reflect reality a bit, we will also approximately define the probabilities of bacteria and the antibiotic results with the prob parameter.
sample_size <- 20000
-data <- data.frame(date = sample(dates, size = sample_size, replace = TRUE),
- patient_id = sample(patients, size = sample_size, replace = TRUE),
- hospital = sample(hospitals, size = sample_size, replace = TRUE,
- prob = c(0.30, 0.35, 0.15, 0.20)),
- bacteria = sample(bacteria, size = sample_size, replace = TRUE,
- prob = c(0.50, 0.25, 0.15, 0.10)),
- AMX = sample(ab_interpretations, size = sample_size, replace = TRUE,
- prob = c(0.60, 0.05, 0.35)),
- AMC = sample(ab_interpretations, size = sample_size, replace = TRUE,
- prob = c(0.75, 0.10, 0.15)),
- CIP = sample(ab_interpretations, size = sample_size, replace = TRUE,
- prob = c(0.80, 0.00, 0.20)),
- GEN = sample(ab_interpretations, size = sample_size, replace = TRUE,
- prob = c(0.92, 0.00, 0.08))
- )sample_size <- 20000
+data <- data.frame(date = sample(dates, size = sample_size, replace = TRUE),
+ patient_id = sample(patients, size = sample_size, replace = TRUE),
+ hospital = sample(hospitals, size = sample_size, replace = TRUE,
+ prob = c(0.30, 0.35, 0.15, 0.20)),
+ bacteria = sample(bacteria, size = sample_size, replace = TRUE,
+ prob = c(0.50, 0.25, 0.15, 0.10)),
+ AMX = sample(ab_interpretations, size = sample_size, replace = TRUE,
+ prob = c(0.60, 0.05, 0.35)),
+ AMC = sample(ab_interpretations, size = sample_size, replace = TRUE,
+ prob = c(0.75, 0.10, 0.15)),
+ CIP = sample(ab_interpretations, size = sample_size, replace = TRUE,
+ prob = c(0.80, 0.00, 0.20)),
+ GEN = sample(ab_interpretations, size = sample_size, replace = TRUE,
+ prob = c(0.92, 0.00, 0.08))
+ )Using the left_join() function from the dplyr package, we can ‘map’ the gender to the patient ID using the patients_table object we created earlier:
The resulting data set contains 20,000 blood culture isolates. With the head() function we can preview the first 6 values of this data set:
| date | @@ -327,70 +327,70 @@||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2011-09-06 | -Z5 | -Hospital B | -Escherichia coli | -R | -S | -S | -S | -F | -||||||||
| 2015-03-21 | -E7 | +2011-07-20 | +X1 | Hospital C | -Escherichia coli | +Staphylococcus aureus | R | I | S | S | -M | -|||||
| 2010-08-11 | -X6 | -Hospital C | -Escherichia coli | -S | -S | -R | -S | F | ||||||||
| 2012-06-16 | -E10 | -Hospital D | -Staphylococcus aureus | -R | +2017-12-18 | +Q8 | +Hospital B | +Streptococcus pneumoniae | S | S | S | -M | +S | +F | ||
| 2016-12-29 | -J3 | -Hospital C | -Escherichia coli | +2015-12-08 | +X5 | +Hospital B | +Staphylococcus aureus | +S | +S | R | -S | -S | -S | -M | +R | +F |
| 2010-04-09 | -Q3 | +2011-12-08 | +X3 | Hospital B | +Klebsiella pneumoniae | +S | +S | +S | +S | +F | +||||||
| 2011-01-02 | +E1 | +Hospital A | Streptococcus pneumoniae | R | S | S | S | -F | +M | +|||||||
| 2012-10-22 | +K1 | +Hospital C | +Streptococcus pneumoniae | +S | +S | +S | +S | +M |
Cleaning the data
Use the frequency table function freq() to look specifically for unique values in any variable. For example, for the gender variable:
# Frequency table of `gender` from `data` (20,000 x 9)
#
# Class: factor (numeric)
@@ -411,82 +411,82 @@
#
# Item Count Percent Cum. Count Cum. Percent
# --- ----- ------- -------- ----------- -------------
-# 1 M 10,408 52.0% 10,408 52.0%
-# 2 F 9,592 48.0% 20,000 100.0%So, we can draw at least two conclusions immediately. From a data scientists perspective, the data looks clean: only values M and F. From a researchers perspective: there are slightly more men. Nothing we didn’t already know.
The data is already quite clean, but we still need to transform some variables. The bacteria column now consists of text, and we want to add more variables based on microbial IDs later on. So, we will transform this column to valid IDs. The mutate() function of the dplyr package makes this really easy:
We also want to transform the antibiotics, because in real life data we don’t know if they are really clean. The as.rsi() function ensures reliability and reproducibility in these kind of variables. The mutate_at() will run the as.rsi() function on defined variables:
Finally, we will apply EUCAST rules on our antimicrobial results. In Europe, most medical microbiological laboratories already apply these rules. Our package features their latest insights on intrinsic resistance and exceptional phenotypes. Moreover, the eucast_rules() function can also apply additional rules, like forcing
Because the amoxicillin (column AMX) and amoxicillin/clavulanic acid (column AMC) in our data were generated randomly, some rows will undoubtedly contain AMX = S and AMC = R, which is technically impossible. The eucast_rules() fixes this:
data <- eucast_rules(data, col_mo = "bacteria")
-#
-# Rules by the European Committee on Antimicrobial Susceptibility Testing (EUCAST)
-# http://eucast.org/
-#
-# EUCAST Clinical Breakpoints (v9.0, 2019)
-# Aerococcus sanguinicola (no new changes)
-# Aerococcus urinae (no new changes)
-# Anaerobic Gram-negatives (no new changes)
-# Anaerobic Gram-positives (no new changes)
-# Campylobacter coli (no new changes)
-# Campylobacter jejuni (no new changes)
-# Enterobacteriales (Order) (no new changes)
-# Enterococcus (no new changes)
-# Haemophilus influenzae (no new changes)
-# Kingella kingae (no new changes)
-# Moraxella catarrhalis (no new changes)
-# Pasteurella multocida (no new changes)
-# Staphylococcus (no new changes)
-# Streptococcus groups A, B, C, G (no new changes)
-# Streptococcus pneumoniae (1,443 new changes)
-# Viridans group streptococci (no new changes)
-#
-# EUCAST Expert Rules, Intrinsic Resistance and Exceptional Phenotypes (v3.1, 2016)
-# Table 01: Intrinsic resistance in Enterobacteriaceae (1,332 new changes)
-# Table 02: Intrinsic resistance in non-fermentative Gram-negative bacteria (no new changes)
-# Table 03: Intrinsic resistance in other Gram-negative bacteria (no new changes)
-# Table 04: Intrinsic resistance in Gram-positive bacteria (2,723 new changes)
-# Table 08: Interpretive rules for B-lactam agents and Gram-positive cocci (no new changes)
-# Table 09: Interpretive rules for B-lactam agents and Gram-negative rods (no new changes)
-# Table 11: Interpretive rules for macrolides, lincosamides, and streptogramins (no new changes)
-# Table 12: Interpretive rules for aminoglycosides (no new changes)
-# Table 13: Interpretive rules for quinolones (no new changes)
-#
-# Other rules
-# Non-EUCAST: amoxicillin/clav acid = S where ampicillin = S (2,213 new changes)
-# Non-EUCAST: ampicillin = R where amoxicillin/clav acid = R (127 new changes)
-# Non-EUCAST: piperacillin = R where piperacillin/tazobactam = R (no new changes)
-# Non-EUCAST: piperacillin/tazobactam = S where piperacillin = S (no new changes)
-# Non-EUCAST: trimethoprim = R where trimethoprim/sulfa = R (no new changes)
-# Non-EUCAST: trimethoprim/sulfa = S where trimethoprim = S (no new changes)
-#
-# --------------------------------------------------------------------------
-# EUCAST rules affected 6,513 out of 20,000 rows, making a total of 7,838 edits
-# => added 0 test results
-#
-# => changed 7,838 test results
-# - 115 test results changed from S to I
-# - 4,719 test results changed from S to R
-# - 1,077 test results changed from I to S
-# - 335 test results changed from I to R
-# - 1,573 test results changed from R to S
-# - 19 test results changed from R to I
-# --------------------------------------------------------------------------
-#
-# Use verbose = TRUE to get a data.frame with all specified edits instead.data <- eucast_rules(data, col_mo = "bacteria")
+#
+# Rules by the European Committee on Antimicrobial Susceptibility Testing (EUCAST)
+# http://eucast.org/
+#
+# EUCAST Clinical Breakpoints (v9.0, 2019)
+# Aerococcus sanguinicola (no new changes)
+# Aerococcus urinae (no new changes)
+# Anaerobic Gram-negatives (no new changes)
+# Anaerobic Gram-positives (no new changes)
+# Campylobacter coli (no new changes)
+# Campylobacter jejuni (no new changes)
+# Enterobacteriales (Order) (no new changes)
+# Enterococcus (no new changes)
+# Haemophilus influenzae (no new changes)
+# Kingella kingae (no new changes)
+# Moraxella catarrhalis (no new changes)
+# Pasteurella multocida (no new changes)
+# Staphylococcus (no new changes)
+# Streptococcus groups A, B, C, G (no new changes)
+# Streptococcus pneumoniae (1,481 new changes)
+# Viridans group streptococci (no new changes)
+#
+# EUCAST Expert Rules, Intrinsic Resistance and Exceptional Phenotypes (v3.1, 2016)
+# Table 01: Intrinsic resistance in Enterobacteriaceae (1,328 new changes)
+# Table 02: Intrinsic resistance in non-fermentative Gram-negative bacteria (no new changes)
+# Table 03: Intrinsic resistance in other Gram-negative bacteria (no new changes)
+# Table 04: Intrinsic resistance in Gram-positive bacteria (2,778 new changes)
+# Table 08: Interpretive rules for B-lactam agents and Gram-positive cocci (no new changes)
+# Table 09: Interpretive rules for B-lactam agents and Gram-negative rods (no new changes)
+# Table 11: Interpretive rules for macrolides, lincosamides, and streptogramins (no new changes)
+# Table 12: Interpretive rules for aminoglycosides (no new changes)
+# Table 13: Interpretive rules for quinolones (no new changes)
+#
+# Other rules
+# Non-EUCAST: amoxicillin/clav acid = S where ampicillin = S (2,245 new changes)
+# Non-EUCAST: ampicillin = R where amoxicillin/clav acid = R (118 new changes)
+# Non-EUCAST: piperacillin = R where piperacillin/tazobactam = R (no new changes)
+# Non-EUCAST: piperacillin/tazobactam = S where piperacillin = S (no new changes)
+# Non-EUCAST: trimethoprim = R where trimethoprim/sulfa = R (no new changes)
+# Non-EUCAST: trimethoprim/sulfa = S where trimethoprim = S (no new changes)
+#
+# --------------------------------------------------------------------------
+# EUCAST rules affected 6,581 out of 20,000 rows, making a total of 7,950 edits
+# => added 0 test results
+#
+# => changed 7,950 test results
+# - 120 test results changed from S to I
+# - 4,812 test results changed from S to R
+# - 1,089 test results changed from I to S
+# - 317 test results changed from I to R
+# - 1,588 test results changed from R to S
+# - 24 test results changed from R to I
+# --------------------------------------------------------------------------
+#
+# Use verbose = TRUE (on your original data) to get a data.frame with all specified edits instead.Adding new variables
Now that we have the microbial ID, we can add some taxonomic properties:
-data <- data %>%
- mutate(gramstain = mo_gramstain(bacteria),
- genus = mo_genus(bacteria),
- species = mo_species(bacteria))data <- data %>%
+ mutate(gramstain = mo_gramstain(bacteria),
+ genus = mo_genus(bacteria),
+ species = mo_species(bacteria))First isolates
@@ -497,23 +497,23 @@(…) When preparing a cumulative antibiogram to guide clinical decisions about empirical antimicrobial therapy of initial infections, only the first isolate of a given species per patient, per analysis period (eg, one year) should be included, irrespective of body site, antimicrobial susceptibility profile, or other phenotypical characteristics (eg, biotype). The first isolate is easily identified, and cumulative antimicrobial susceptibility test data prepared using the first isolate are generally comparable to cumulative antimicrobial susceptibility test data calculated by other methods, providing duplicate isolates are excluded.
M39-A4 Analysis and Presentation of Cumulative Antimicrobial Susceptibility Test Data, 4th Edition. CLSI, 2014. Chapter 6.4
This AMR package includes this methodology with the first_isolate() function. It adopts the episode of a year (can be changed by user) and it starts counting days after every selected isolate. This new variable can easily be added to our data:
data <- data %>%
- mutate(first = first_isolate(.))
-# NOTE: Using column `bacteria` as input for `col_mo`.
-# NOTE: Using column `date` as input for `col_date`.
-# NOTE: Using column `patient_id` as input for `col_patient_id`.
-# => Found 5,719 first isolates (28.6% of total)So only 28.6% is suitable for resistance analysis! We can now filter on it with the filter() function, also from the dplyr package:
data <- data %>%
+ mutate(first = first_isolate(.))
+# NOTE: Using column `bacteria` as input for `col_mo`.
+# NOTE: Using column `date` as input for `col_date`.
+# NOTE: Using column `patient_id` as input for `col_patient_id`.
+# => Found 5,679 first isolates (28.4% of total)So only 28.4% is suitable for resistance analysis! We can now filter on it with the filter() function, also from the dplyr package:
For future use, the above two syntaxes can be shortened with the filter_first_isolate() function:
First weighted isolates
-We made a slight twist to the CLSI algorithm, to take into account the antimicrobial susceptibility profile. Have a look at all isolates of patient S7, sorted on date:
+We made a slight twist to the CLSI algorithm, to take into account the antimicrobial susceptibility profile. Have a look at all isolates of patient B6, sorted on date:
| isolate | @@ -529,30 +529,30 @@|||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -2010-01-28 | -S7 | +2010-02-26 | +B6 | B_ESCHR_COL | R | -I | S | +R | S | TRUE | ||
| 2 | -2010-02-07 | -S7 | +2010-05-20 | +B6 | B_ESCHR_COL | -S | -S | -S | R | +S | +S | +S | FALSE |
| 3 | -2010-03-16 | -S7 | +2010-05-28 | +B6 | B_ESCHR_COL | R | S | @@ -562,10 +562,10 @@||||||
| 4 | -2010-10-09 | -S7 | +2010-06-27 | +B6 | B_ESCHR_COL | -S | +R | S | S | S | @@ -573,8 +573,8 @@|||
| 5 | -2011-01-25 | -S7 | +2010-09-06 | +B6 | B_ESCHR_COL | R | S | @@ -584,19 +584,19 @@||||||
| 6 | -2011-02-16 | -S7 | +2010-10-16 | +B6 | B_ESCHR_COL | +R | S | S | S | -S | -TRUE | +FALSE | |
| 7 | -2011-02-24 | -S7 | +2010-10-20 | +B6 | B_ESCHR_COL | S | S | @@ -606,32 +606,32 @@||||||
| 8 | -2011-03-30 | -S7 | +2010-11-01 | +B6 | B_ESCHR_COL | -R | +S | +S | S | R | -S | FALSE | |
| 9 | -2011-04-25 | -S7 | +2010-11-14 | +B6 | B_ESCHR_COL | S | S | S | -R | +S | FALSE | ||
| 10 | -2011-05-06 | -S7 | +2011-01-25 | +B6 | B_ESCHR_COL | -S | +R | S | S | S | @@ -639,18 +639,18 @@
Only 2 isolates are marked as ‘first’ according to CLSI guideline. But when reviewing the antibiogram, it is obvious that some isolates are absolutely different strains and should be included too. This is why we weigh isolates, based on their antibiogram. The key_antibiotics() function adds a vector with 18 key antibiotics: 6 broad spectrum ones, 6 small spectrum for Gram negatives and 6 small spectrum for Gram positives. These can be defined by the user.
Only 1 isolates are marked as ‘first’ according to CLSI guideline. But when reviewing the antibiogram, it is obvious that some isolates are absolutely different strains and should be included too. This is why we weigh isolates, based on their antibiogram. The key_antibiotics() function adds a vector with 18 key antibiotics: 6 broad spectrum ones, 6 small spectrum for Gram negatives and 6 small spectrum for Gram positives. These can be defined by the user.
If a column exists with a name like ‘key(…)ab’ the first_isolate() function will automatically use it and determine the first weighted isolates. Mind the NOTEs in below output:
data <- data %>%
- mutate(keyab = key_antibiotics(.)) %>%
- mutate(first_weighted = first_isolate(.))
-# NOTE: Using column `bacteria` as input for `col_mo`.
-# NOTE: Using column `bacteria` as input for `col_mo`.
-# NOTE: Using column `date` as input for `col_date`.
-# NOTE: Using column `patient_id` as input for `col_patient_id`.
-# NOTE: Using column `keyab` as input for `col_keyantibiotics`. Use col_keyantibiotics = FALSE to prevent this.
-# [Criterion] Inclusion based on key antibiotics, ignoring I.
-# => Found 15,097 first weighted isolates (75.5% of total)data <- data %>%
+ mutate(keyab = key_antibiotics(.)) %>%
+ mutate(first_weighted = first_isolate(.))
+# NOTE: Using column `bacteria` as input for `col_mo`.
+# NOTE: Using column `bacteria` as input for `col_mo`.
+# NOTE: Using column `date` as input for `col_date`.
+# NOTE: Using column `patient_id` as input for `col_patient_id`.
+# NOTE: Using column `keyab` as input for `col_keyantibiotics`. Use col_keyantibiotics = FALSE to prevent this.
+# [Criterion] Inclusion based on key antibiotics, ignoring I.
+# => Found 15,003 first weighted isolates (75.0% of total)| isolate | @@ -667,118 +667,118 @@||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -2010-01-28 | -S7 | +2010-02-26 | +B6 | B_ESCHR_COL | R | -I | S | +R | S | TRUE | TRUE | ||
| 2 | -2010-02-07 | -S7 | +2010-05-20 | +B6 | B_ESCHR_COL | -S | -S | -S | R | +S | +S | +S | FALSE | TRUE |
| 3 | -2010-03-16 | -S7 | +2010-05-28 | +B6 | B_ESCHR_COL | R | S | S | S | FALSE | -TRUE | +FALSE | ||
| 4 | -2010-10-09 | -S7 | +2010-06-27 | +B6 | B_ESCHR_COL | -S | +R | S | S | S | FALSE | -TRUE | +FALSE | |
| 5 | -2011-01-25 | -S7 | +2010-09-06 | +B6 | B_ESCHR_COL | R | S | S | S | FALSE | -TRUE | +FALSE | ||
| 6 | -2011-02-16 | -S7 | +2010-10-16 | +B6 | B_ESCHR_COL | +R | S | S | S | -S | -TRUE | -TRUE | +FALSE | +FALSE |
| 7 | -2011-02-24 | -S7 | +2010-10-20 | +B6 | B_ESCHR_COL | S | S | S | S | FALSE | -FALSE | +TRUE | ||
| 8 | -2011-03-30 | -S7 | +2010-11-01 | +B6 | B_ESCHR_COL | -R | +S | +S | S | R | -S | FALSE | TRUE | |
| 9 | -2011-04-25 | -S7 | +2010-11-14 | +B6 | B_ESCHR_COL | S | S | S | -R | +S | FALSE | TRUE | ||
| 10 | -2011-05-06 | -S7 | +2011-01-25 | +B6 | B_ESCHR_COL | -S | +R | S | S | S | @@ -787,18 +787,19 @@
Instead of 2, now 9 isolates are flagged. In total, 75.5% of all isolates are marked ‘first weighted’ - 46.9% more than when using the CLSI guideline. In real life, this novel algorithm will yield 5-10% more isolates than the classic CLSI guideline.
+Instead of 1, now 6 isolates are flagged. In total, 75% of all isolates are marked ‘first weighted’ - 46.6% more than when using the CLSI guideline. In real life, this novel algorithm will yield 5-10% more isolates than the classic CLSI guideline.
As with filter_first_isolate(), there’s a shortcut for this new algorithm too:
So we end up with 15,097 isolates for analysis.
+ +So we end up with 15,003 isolates for analysis.
We can remove unneeded columns:
- +Now our data looks like:
- +
(omitted 29 entries, n = 57 [11.4%])
-
-# our transformed antibiotic columns
-# amoxicillin/clavulanic acid (J01CR02) as an example
-data %>% freq(AMC_ND2)
+# our transformed antibiotic columns
+# amoxicillin/clavulanic acid (J01CR02) as an example
+data %>% freq(AMC_ND2)Frequency table of AMC_ND2 from data (500 x 54)
Class: factor > ordered > rsi (numeric)
Length: 500 (of which NA: 19 = 3.80%)
diff --git a/docs/articles/benchmarks.html b/docs/articles/benchmarks.html
index d950825d..84caae51 100644
--- a/docs/articles/benchmarks.html
+++ b/docs/articles/benchmarks.html
@@ -40,7 +40,7 @@
Benchmarks
Matthijs S. Berends
-01 July 2019
+10 July 2019
benchmarks.RmdOne of the most important features of this package is the complete microbial taxonomic database, supplied by the Catalogue of Life. We created a function as.mo() that transforms any user input value to a valid microbial ID by using intelligent rules combined with the taxonomic tree of Catalogue of Life.
Using the microbenchmark package, we can review the calculation performance of this function. Its function microbenchmark() runs different input expressions independently of each other and measures their time-to-result.
In the next test, we try to ‘coerce’ different input values for Staphylococcus aureus. The actual result is the same every time: it returns its MO code B_STPHY_AUR (B stands for Bacteria, the taxonomic kingdom).
But the calculation time differs a lot:
-S.aureus <- microbenchmark(as.mo("sau"),
- as.mo("stau"),
- as.mo("staaur"),
- as.mo("STAAUR"),
- as.mo("S. aureus"),
- as.mo("S. aureus"),
- as.mo("Staphylococcus aureus"),
- times = 10)
-print(S.aureus, unit = "ms", signif = 2)
-# Unit: milliseconds
-# expr min lq mean median uq max neval
-# as.mo("sau") 18.0 18.0 22 18.0 18.0 61 10
-# as.mo("stau") 65.0 65.0 70 66.0 66.0 110 10
-# as.mo("staaur") 18.0 18.0 33 18.0 62.0 81 10
-# as.mo("STAAUR") 18.0 18.0 18 18.0 18.0 19 10
-# as.mo("S. aureus") 52.0 52.0 61 52.0 53.0 97 10
-# as.mo("S. aureus") 52.0 52.0 71 53.0 97.0 150 10
-# as.mo("Staphylococcus aureus") 8.1 8.1 14 8.1 8.2 63 10S.aureus <- microbenchmark(as.mo("sau"),
+ as.mo("stau"),
+ as.mo("staaur"),
+ as.mo("STAAUR"),
+ as.mo("S. aureus"),
+ as.mo("S. aureus"),
+ as.mo("Staphylococcus aureus"),
+ times = 10)
+print(S.aureus, unit = "ms", signif = 2)
+# Unit: milliseconds
+# expr min lq mean median uq max neval
+# as.mo("sau") 8.5 8.7 12.0 8.9 9.4 26 10
+# as.mo("stau") 31.0 32.0 42.0 33.0 34.0 120 10
+# as.mo("staaur") 8.6 8.7 11.0 9.1 9.2 26 10
+# as.mo("STAAUR") 8.7 9.1 9.3 9.2 9.4 11 10
+# as.mo("S. aureus") 23.0 23.0 30.0 24.0 40.0 46 10
+# as.mo("S. aureus") 22.0 23.0 27.0 24.0 25.0 41 10
+# as.mo("Staphylococcus aureus") 3.9 4.0 5.7 4.1 4.4 20 10In the table above, all measurements are shown in milliseconds (thousands of seconds). A value of 5 milliseconds means it can determine 200 input values per second. It case of 100 milliseconds, this is only 10 input values per second. The second input is the only one that has to be looked up thoroughly. All the others are known codes (the first one is a WHONET code) or common laboratory codes, or common full organism names like the last one. Full organism names are always preferred.
To achieve this speed, the as.mo function also takes into account the prevalence of human pathogenic microorganisms. The downside is of course that less prevalent microorganisms will be determined less fast. See this example for the ID of Thermus islandicus (B_THERMS_ISL), a bug probably never found before in humans:
T.islandicus <- microbenchmark(as.mo("theisl"),
- as.mo("THEISL"),
- as.mo("T. islandicus"),
- as.mo("T. islandicus"),
- as.mo("Thermus islandicus"),
- times = 10)
-print(T.islandicus, unit = "ms", signif = 2)
-# Unit: milliseconds
-# expr min lq mean median uq max neval
-# as.mo("theisl") 390 390 420 440 440 440 10
-# as.mo("THEISL") 390 390 410 400 440 440 10
-# as.mo("T. islandicus") 210 210 230 220 250 270 10
-# as.mo("T. islandicus") 210 210 240 260 260 280 10
-# as.mo("Thermus islandicus") 72 72 92 73 120 130 10That takes 6.8 times as much time on average. A value of 100 milliseconds means it can only determine ~10 different input values per second. We can conclude that looking up arbitrary codes of less prevalent microorganisms is the worst way to go, in terms of calculation performance. Full names (like Thermus islandicus) are almost fast - these are the most probable input from most data sets.
+T.islandicus <- microbenchmark(as.mo("theisl"),
+ as.mo("THEISL"),
+ as.mo("T. islandicus"),
+ as.mo("T. islandicus"),
+ as.mo("Thermus islandicus"),
+ times = 10)
+print(T.islandicus, unit = "ms", signif = 2)
+# Unit: milliseconds
+# expr min lq mean median uq max neval
+# as.mo("theisl") 290 310 320 320 320 350 10
+# as.mo("THEISL") 290 300 310 310 330 340 10
+# as.mo("T. islandicus") 140 140 150 150 160 170 10
+# as.mo("T. islandicus") 140 140 150 150 170 180 10
+# as.mo("Thermus islandicus") 50 52 58 54 68 70 10That takes 10.2 times as much time on average. A value of 100 milliseconds means it can only determine ~10 different input values per second. We can conclude that looking up arbitrary codes of less prevalent microorganisms is the worst way to go, in terms of calculation performance. Full names (like Thermus islandicus) are almost fast - these are the most probable input from most data sets.
In the figure below, we compare Escherichia coli (which is very common) with Prevotella brevis (which is moderately common) and with Thermus islandicus (which is very uncommon):
-par(mar = c(5, 16, 4, 2)) # set more space for left margin text (16)
-
-boxplot(microbenchmark(as.mo("Thermus islandicus"),
- as.mo("Prevotella brevis"),
- as.mo("Escherichia coli"),
- as.mo("T. islandicus"),
- as.mo("P. brevis"),
- as.mo("E. coli"),
- times = 10),
- horizontal = TRUE, las = 1, unit = "s", log = FALSE,
- xlab = "", ylab = "Time in seconds",
- main = "Benchmarks per prevalence")par(mar = c(5, 16, 4, 2)) # set more space for left margin text (16)
+
+boxplot(microbenchmark(as.mo("Thermus islandicus"),
+ as.mo("Prevotella brevis"),
+ as.mo("Escherichia coli"),
+ as.mo("T. islandicus"),
+ as.mo("P. brevis"),
+ as.mo("E. coli"),
+ times = 10),
+ horizontal = TRUE, las = 1, unit = "s", log = FALSE,
+ xlab = "", ylab = "Time in seconds",
+ main = "Benchmarks per prevalence")
Uncommon microorganisms take a lot more time than common microorganisms. To relieve this pitfall and further improve performance, two important calculations take almost no time at all: repetitive results and already precalculated results.
Repetitive results
Repetitive results are unique values that are present more than once. Unique values will only be calculated once by as.mo(). We will use mo_fullname() for this test - a helper function that returns the full microbial name (genus, species and possibly subspecies) which uses as.mo() internally.
library(dplyr)
-# take all MO codes from the septic_patients data set
-x <- septic_patients$mo %>%
- # keep only the unique ones
- unique() %>%
- # pick 50 of them at random
- sample(50) %>%
- # paste that 10,000 times
- rep(10000) %>%
- # scramble it
- sample()
-
-# got indeed 50 times 10,000 = half a million?
-length(x)
-# [1] 500000
-
-# and how many unique values do we have?
-n_distinct(x)
-# [1] 50
-
-# now let's see:
-run_it <- microbenchmark(mo_fullname(x),
- times = 10)
-print(run_it, unit = "ms", signif = 3)
-# Unit: milliseconds
-# expr min lq mean median uq max neval
-# mo_fullname(x) 1050 1050 1100 1090 1120 1230 10So transforming 500,000 values (!!) of 50 unique values only takes 1.09 seconds (1092 ms). You only lose time on your unique input values.
+library(dplyr)
+# take all MO codes from the septic_patients data set
+x <- septic_patients$mo %>%
+ # keep only the unique ones
+ unique() %>%
+ # pick 50 of them at random
+ sample(50) %>%
+ # paste that 10,000 times
+ rep(10000) %>%
+ # scramble it
+ sample()
+
+# got indeed 50 times 10,000 = half a million?
+length(x)
+# [1] 500000
+
+# and how many unique values do we have?
+n_distinct(x)
+# [1] 50
+
+# now let's see:
+run_it <- microbenchmark(mo_fullname(x),
+ times = 10)
+print(run_it, unit = "ms", signif = 3)
+# Unit: milliseconds
+# expr min lq mean median uq max neval
+# mo_fullname(x) 611 628 643 635 650 714 10So transforming 500,000 values (!!) of 50 unique values only takes 0.63 seconds (634 ms). You only lose time on your unique input values.
Precalculated results
What about precalculated results? If the input is an already precalculated result of a helper function like mo_fullname(), it almost doesn’t take any time at all (see ‘C’ below):
run_it <- microbenchmark(A = mo_fullname("B_STPHY_AUR"),
- B = mo_fullname("S. aureus"),
- C = mo_fullname("Staphylococcus aureus"),
- times = 10)
-print(run_it, unit = "ms", signif = 3)
-# Unit: milliseconds
-# expr min lq mean median uq max neval
-# A 13.00 13.20 13.60 13.60 14.00 14.40 10
-# B 49.40 50.00 57.50 51.90 52.40 103.00 10
-# C 1.52 1.72 1.81 1.78 1.98 1.99 10So going from mo_fullname("Staphylococcus aureus") to "Staphylococcus aureus" takes 0.0018 seconds - it doesn’t even start calculating if the result would be the same as the expected resulting value. That goes for all helper functions:
run_it <- microbenchmark(A = mo_species("aureus"),
- B = mo_genus("Staphylococcus"),
- C = mo_fullname("Staphylococcus aureus"),
- D = mo_family("Staphylococcaceae"),
- E = mo_order("Bacillales"),
- F = mo_class("Bacilli"),
- G = mo_phylum("Firmicutes"),
- H = mo_kingdom("Bacteria"),
- times = 10)
-print(run_it, unit = "ms", signif = 3)
-# Unit: milliseconds
-# expr min lq mean median uq max neval
-# A 0.612 0.623 0.685 0.653 0.789 0.814 10
-# B 0.556 0.575 0.680 0.671 0.689 0.958 10
-# C 1.520 1.710 1.800 1.820 1.950 1.970 10
-# D 0.547 0.665 0.723 0.688 0.811 0.997 10
-# E 0.490 0.541 0.633 0.629 0.748 0.756 10
-# F 0.482 0.569 0.612 0.590 0.663 0.756 10
-# G 0.551 0.558 0.601 0.586 0.632 0.735 10
-# H 0.494 0.564 0.595 0.575 0.608 0.757 10run_it <- microbenchmark(A = mo_fullname("B_STPHY_AUR"),
+ B = mo_fullname("S. aureus"),
+ C = mo_fullname("Staphylococcus aureus"),
+ times = 10)
+print(run_it, unit = "ms", signif = 3)
+# Unit: milliseconds
+# expr min lq mean median uq max neval
+# A 6.730 7.030 8.030 7.750 8.72 9.73 10
+# B 22.400 23.000 27.100 23.600 27.10 46.00 10
+# C 0.835 0.877 0.978 0.925 1.12 1.18 10So going from mo_fullname("Staphylococcus aureus") to "Staphylococcus aureus" takes 0.0009 seconds - it doesn’t even start calculating if the result would be the same as the expected resulting value. That goes for all helper functions:
run_it <- microbenchmark(A = mo_species("aureus"),
+ B = mo_genus("Staphylococcus"),
+ C = mo_fullname("Staphylococcus aureus"),
+ D = mo_family("Staphylococcaceae"),
+ E = mo_order("Bacillales"),
+ F = mo_class("Bacilli"),
+ G = mo_phylum("Firmicutes"),
+ H = mo_kingdom("Bacteria"),
+ times = 10)
+print(run_it, unit = "ms", signif = 3)
+# Unit: milliseconds
+# expr min lq mean median uq max neval
+# A 0.468 0.470 0.533 0.489 0.595 0.690 10
+# B 0.504 0.513 0.555 0.520 0.571 0.711 10
+# C 0.629 0.687 0.864 0.855 1.050 1.130 10
+# D 0.505 0.515 0.575 0.530 0.649 0.767 10
+# E 0.442 0.457 0.529 0.481 0.531 0.774 10
+# F 0.447 0.510 0.554 0.568 0.609 0.618 10
+# G 0.443 0.470 0.492 0.477 0.506 0.601 10
+# H 0.448 0.459 0.491 0.466 0.515 0.633 10Of course, when running mo_phylum("Firmicutes") the function has zero knowledge about the actual microorganism, namely S. aureus. But since the result would be "Firmicutes" too, there is no point in calculating the result. And because this package ‘knows’ all phyla of all known bacteria (according to the Catalogue of Life), it can just return the initial value immediately.
Results in other languages
When the system language is non-English and supported by this AMR package, some functions will have a translated result. This almost does’t take extra time:
mo_fullname("CoNS", language = "en") # or just mo_fullname("CoNS") on an English system
-# [1] "Coagulase-negative Staphylococcus (CoNS)"
-
-mo_fullname("CoNS", language = "es") # or just mo_fullname("CoNS") on a Spanish system
-# [1] "Staphylococcus coagulasa negativo (SCN)"
-
-mo_fullname("CoNS", language = "nl") # or just mo_fullname("CoNS") on a Dutch system
-# [1] "Coagulase-negatieve Staphylococcus (CNS)"
-
-run_it <- microbenchmark(en = mo_fullname("CoNS", language = "en"),
- de = mo_fullname("CoNS", language = "de"),
- nl = mo_fullname("CoNS", language = "nl"),
- es = mo_fullname("CoNS", language = "es"),
- it = mo_fullname("CoNS", language = "it"),
- fr = mo_fullname("CoNS", language = "fr"),
- pt = mo_fullname("CoNS", language = "pt"),
- times = 10)
-print(run_it, unit = "ms", signif = 4)
-# Unit: milliseconds
-# expr min lq mean median uq max neval
-# en 43.00 43.12 45.51 44.82 44.89 56.61 10
-# de 46.47 46.99 52.11 47.57 48.11 93.77 10
-# nl 60.86 62.72 67.57 63.69 63.99 108.20 10
-# es 45.74 46.05 52.37 46.42 47.98 103.00 10
-# it 45.84 45.89 51.90 47.66 47.73 94.83 10
-# fr 45.97 46.92 47.44 47.76 47.86 48.49 10
-# pt 45.93 46.77 47.36 47.77 47.93 48.12 10mo_fullname("CoNS", language = "en") # or just mo_fullname("CoNS") on an English system
+# [1] "Coagulase-negative Staphylococcus (CoNS)"
+
+mo_fullname("CoNS", language = "es") # or just mo_fullname("CoNS") on a Spanish system
+# [1] "Staphylococcus coagulasa negativo (SCN)"
+
+mo_fullname("CoNS", language = "nl") # or just mo_fullname("CoNS") on a Dutch system
+# [1] "Coagulase-negatieve Staphylococcus (CNS)"
+
+run_it <- microbenchmark(en = mo_fullname("CoNS", language = "en"),
+ de = mo_fullname("CoNS", language = "de"),
+ nl = mo_fullname("CoNS", language = "nl"),
+ es = mo_fullname("CoNS", language = "es"),
+ it = mo_fullname("CoNS", language = "it"),
+ fr = mo_fullname("CoNS", language = "fr"),
+ pt = mo_fullname("CoNS", language = "pt"),
+ times = 10)
+print(run_it, unit = "ms", signif = 4)
+# Unit: milliseconds
+# expr min lq mean median uq max neval
+# en 19.29 20.39 22.42 20.64 21.26 38.98 10
+# de 22.03 22.40 23.78 23.08 23.60 31.74 10
+# nl 27.98 28.60 29.22 28.87 30.13 30.53 10
+# es 21.44 22.97 30.04 24.11 45.71 46.46 10
+# it 21.17 21.82 22.47 22.44 23.27 23.52 10
+# fr 20.75 21.58 24.10 22.04 22.41 42.96 10
+# pt 21.24 21.92 24.31 22.75 23.26 39.91 10Currently supported are German, Dutch, Spanish, Italian, French and Portuguese.
How to create frequency tables
Matthijs S. Berends
-01 July 2019
+10 July 2019
freq.RmdFrequencies of one variable
To only show and quickly review the content of one variable, you can just select this variable in various ways. Let’s say we want to get the frequencies of the gender variable of the septic_patients dataset:
# Any of these will work:
-# freq(septic_patients$gender)
-# freq(septic_patients[, "gender"])
-
-# Using tidyverse:
-# septic_patients$gender %>% freq()
-# septic_patients[, "gender"] %>% freq()
-# septic_patients %>% freq("gender")
-
-# Probably the fastest and easiest:
-septic_patients %>% freq(gender) # Any of these will work:
+# freq(septic_patients$gender)
+# freq(septic_patients[, "gender"])
+
+# Using tidyverse:
+# septic_patients$gender %>% freq()
+# septic_patients[, "gender"] %>% freq()
+# septic_patients %>% freq("gender")
+
+# Probably the fastest and easiest:
+septic_patients %>% freq(gender) Frequency table of gender from septic_patients (2,000 x 49)
Class: character
Length: 2,000 (of which NA: 0 = 0.00%)
@@ -262,22 +262,22 @@ Longest: 1
Multiple variables will be pasted into one variable to review individual cases, keeping a univariate frequency table.
For illustration, we could add some more variables to the septic_patients dataset to learn about bacterial properties:
Now all variables of the microorganisms dataset have been joined to the septic_patients dataset. The microorganisms dataset consists of the following variables:
colnames(microorganisms)
-# [1] "mo" "col_id" "fullname" "kingdom" "phylum"
-# [6] "class" "order" "family" "genus" "species"
-# [11] "subspecies" "rank" "ref" "species_id" "source"
-# [16] "prevalence"colnames(microorganisms)
+# [1] "mo" "col_id" "fullname" "kingdom" "phylum"
+# [6] "class" "order" "family" "genus" "species"
+# [11] "subspecies" "rank" "ref" "species_id" "source"
+# [16] "prevalence"If we compare the dimensions between the old and new dataset, we can see that these 15 variables were added:
- +So now the genus and species variables are available. A frequency table of these combined variables can be created like this:
Frequency table of genus and species from my_patients (2,000 x 64)
Columns: 2
Length: 2,000 (of which NA: 0 = 0.00%)
@@ -423,10 +423,10 @@ Longest: 34
Frequency tables can be created of any input.
In case of numeric values (like integers, doubles, etc.) additional descriptive statistics will be calculated and shown into the header:
-# # get age distribution of unique patients
-septic_patients %>%
- distinct(patient_id, .keep_all = TRUE) %>%
- freq(age, nmax = 5, header = TRUE)# # get age distribution of unique patients
+septic_patients %>%
+ distinct(patient_id, .keep_all = TRUE) %>%
+ freq(age, nmax = 5, header = TRUE)Frequency table of age from a data.frame (981 x 49)
Class: numeric
Length: 981 (of which NA: 0 = 0.00%)
@@ -506,8 +506,8 @@ Outliers: 15 (unique count: 12)
To sort frequencies of factors on their levels instead of item count, use the sort.count parameter.
sort.count is TRUE by default. Compare this default behaviour…
Frequency table of hospital_id from septic_patients (2,000 x 49)
Class: factor (numeric)
Length: 2,000 (of which NA: 0 = 0.00%)
@@ -558,8 +558,8 @@ Unique: 4
… to this, where items are now sorted on factor levels:
- +Frequency table of hospital_id from septic_patients (2,000 x 49)
Class: factor (numeric)
Length: 2,000 (of which NA: 0 = 0.00%)
@@ -610,8 +610,8 @@ Unique: 4
All classes will be printed into the header. Variables with the new rsi class of this AMR package are actually ordered factors and have three classes (look at Class in the header):
Frequency table of AMX from septic_patients (2,000 x 49)
Class: factor > ordered > rsi (numeric)
Length: 2,000 (of which NA: 771 = 38.55%)
@@ -661,8 +661,8 @@ Group: Beta-lactams/penicillins
Frequencies of dates
Frequencies of dates will show the oldest and newest date in the data, and the amount of days between them:
- +Frequency table of date from septic_patients (2,000 x 49)
Class: Date (numeric)
Length: 2,000 (of which NA: 0 = 0.00%)
@@ -728,11 +728,11 @@ Median: 31 July 2009 (47.39%)
Assigning a frequency table to an object
A frequency table is actually a regular data.frame, with the exception that it contains an additional class.
[1] “freq” “data.frame”
Because of this additional class, a frequency table prints like the examples above. But the object itself contains the complete table without a row limitation:
- +[1] 74 5
+
