mirror of
https://github.com/msberends/AMR.git
synced 2026-11-16 00:10:57 +01:00
(v3.1.0.9010) Fix the systematic underestimation of WISCA coverage by only including pathogens that occur in a group and were tested for the regimen, and warn when the prior weight of rare pathogens is large (fixes #337)
This commit is contained in:
dr. M.S. (Matthijs) Berends committed
1 parent
92bb8931c5
commit
fdfdadae14
6 files changed
+93
-13
No files matched your search
+1
-1
@@ -1,5 +1,5 @@
|
||||
Package: AMR
|
||||
Version: 3.1.0.9009
|
||||
Version: 3.1.0.9010
|
||||
Date: 2026-10-30
|
||||
Title: Antimicrobial Resistance Data Analysis
|
||||
Description: Functions to simplify and standardise antimicrobial resistance (AMR)
|
||||
|
||||
@@ -1,6 +1,12 @@
|
||||
# AMR 3.1.0.9009
|
||||
# AMR 3.1.0.9010
|
||||
|
||||
This version fixes a critical bug in WISCA that biased all coverage estimates of `wisca()` and `antibiogram(..., wisca = TRUE)` downwards since AMR 3.0.0. Users are strongly advised to re-run their WISCA analyses.
|
||||
|
||||
### Changed
|
||||
* WISCA (`wisca()` and `antibiogram(..., wisca = TRUE)`)
|
||||
* Fixed a critical bug that systematically underestimated the coverage of every WISCA since AMR 3.0.0. Pathogens absent from a syndromic group, and pathogens without any test result for a regimen, were still included in the model with only prior information: each added a prior weight of 1 to the pathogen incidence, and untested pathogens kept their full incidence weight while counting as about 50% susceptible. The bias was largest for stratified analyses (e.g., up to 14 percentage points for meropenem in random strata of `example_isolates`). Now, only pathogens that occur in a group and have at least one test result for the regimen are included. Please re-run all WISCA analyses (#337)
|
||||
* `retrieve_wisca_parameters()` now returns only the parameters actually used in the model, and `n_total` in the `long_numeric` attribute now counts only isolates of the pathogens included for that regimen
|
||||
* A new warning is given when the prior weight of 1 per pathogen (Bielicki *et al.*, 2016) accounts for more than 5% of the pathogen incidence, as many rare pathogens can then strongly lower the estimated coverage; consider `top_n_microorganisms()` in that case
|
||||
* Microorganism coercion (`as.mo()`)
|
||||
* Fixed the misaligned arrows in the message on taxonomically renamed microorganisms in colour-supporting consoles, for old names whose author reference lacks "*et al.*" (#336)
|
||||
|
||||
|
||||
+24
-2
@@ -156,6 +156,10 @@
|
||||
#'
|
||||
#' **Prior distributions:** Pathogen incidence uses a non-informative \eqn{Dirichlet(1, 1, \ldots, 1)} prior. Susceptibility proportions use the Jeffreys prior, \eqn{\beta(0.5, 0.5)}, except for pathogen-drug combinations with known intrinsic resistance, which use a strongly informative \eqn{\beta(1, 9999)} prior that forces near-zero susceptibility regardless of observed data. Intrinsic resistance is determined using the [intrinsic_resistant] data set, which is based on `r format_eucast_version_nr(names(EUCAST_VERSION_EXPECTED_PHENOTYPES[1]))`.
|
||||
#'
|
||||
#' **Pathogens included in the model:** For each regimen (and group), only pathogens that occur in the data and have at least one test result for the regimen are included. Pathogens without any test result for a regimen do not contribute to its coverage, nor to its `n_total` in `long_numeric`. Use [retrieve_wisca_parameters()] to get the parameters that were used.
|
||||
#'
|
||||
#' **Rare pathogens:** Since every pathogen receives a prior weight of 1 in the Dirichlet distribution, pathogens found only once or a few times are given relatively much weight, which pulls the incidence towards a uniform distribution. Their very uncertain susceptibility can then strongly lower the estimated coverage, without reflecting the clinical expectation. A warning is therefore given when this prior weight accounts for more than 5% of the pathogen incidence in any regimen (or group). Consider limiting the analysis to the most common pathogens using [top_n_microorganisms()], or grouping rare pathogens before running WISCA.
|
||||
#'
|
||||
#' **Interpreting the output:** Overlapping credible intervals between regimens indicate no significant difference in coverage; if a narrower-spectrum regimen overlaps with a broader one, the narrower-spectrum option may be preferred on stewardship grounds. Non-overlapping intervals indicate a clinically meaningful difference. For small sample sizes, consider pooling data from multiple sites to improve precision, provided pathogen distributions are sufficiently similar (Bielicki *et al.*, 2016).
|
||||
#'
|
||||
#' For the full mathematical derivation and worked examples, see the [WISCA vignette](https://amr-for-r.org/articles/WISCA.html).
|
||||
@@ -744,7 +748,8 @@ antibiogram.default <- function(x,
|
||||
coverage = NA_real_,
|
||||
lower_ci = NA_real_,
|
||||
upper_ci = NA_real_,
|
||||
n_total = sum(n_total, na.rm = TRUE),
|
||||
# only pathogens with test results for the regimen enter the WISCA model, see below
|
||||
n_total = sum(n_total[which(n_tested > 0)]),
|
||||
n_tested = sum(n_tested, na.rm = TRUE),
|
||||
n_susceptible = sum(n_susceptible, na.rm = TRUE)
|
||||
)
|
||||
@@ -761,7 +766,24 @@ antibiogram.default <- function(x,
|
||||
out_wisca$group <- out_wisca$ab
|
||||
}
|
||||
|
||||
wisca_parameters <- out
|
||||
# only pathogens that occur in the group (n_total > 0) and have at least one test result for the regimen (n_tested > 0) may contribute to WISCA; otherwise, absent or untested pathogens would enter with prior-only parameters and bias the coverage downwards (#337)
|
||||
wisca_parameters <- out[which(out$n_total > 0 & out$n_tested > 0), , drop = FALSE]
|
||||
|
||||
# the Dirichlet prior of 1 per pathogen (Bielicki et al., 2016) weighs heavily when many pathogens are rare
|
||||
prior_share <- vapply(
|
||||
split(wisca_parameters$n_total, as.character(wisca_parameters$group)),
|
||||
function(n) length(n) / (length(n) + sum(n)),
|
||||
FUN.VALUE = numeric(1)
|
||||
)
|
||||
if (any(prior_share > 0.05)) {
|
||||
warning_(
|
||||
"In ", sum(prior_share > 0.05), " of ", length(prior_share), " WISCA regimen(s) or group(s), the prior weight of 1 per pathogen accounts for more than 5% of the pathogen incidence (up to ", round(max(prior_share) * 100, 1), "%), since many pathogens occur only rarely. ",
|
||||
"Because WISCA gives every pathogen a prior weight of 1 (Bielicki et al., 2016), such rare pathogens can strongly lower the estimated coverage without reflecting the clinical expectation. ",
|
||||
"Consider limiting the analysis to the most common pathogens with {.help [{.fun top_n_microorganisms}](AMR::top_n_microorganisms)}, or grouping rare pathogens before running WISCA.",
|
||||
call = FALSE
|
||||
)
|
||||
}
|
||||
|
||||
wisca_draws <- list()
|
||||
wisca_components <- list()
|
||||
|
||||
|
||||
@@ -59,8 +59,8 @@ times](https://scholar.google.com/citations?view_op=view_citation&hl=en&citation
|
||||
in scientific research.
|
||||
|
||||
After installing this package, R knows [**~103 000 distinct microbial
|
||||
species**](./reference/microorganisms.html) (updated April 2026) and all
|
||||
[**~630 antimicrobial and antiviral
|
||||
species**](./reference/microorganisms.html) (updated October 2026) and
|
||||
all [**~630 antimicrobial and antiviral
|
||||
drugs**](./reference/antimicrobials.html) by name and code (including
|
||||
ATC, EARS-Net, ASIARS-Net, PubChem, LOINC and SNOMED CT), and knows all
|
||||
about valid SIR and MIC values. The integral clinical breakpoint
|
||||
@@ -224,7 +224,7 @@ wisca(example_isolates,
|
||||
|
||||
| Piperacillin/tazobactam | Piperacillin/tazobactam + Gentamicin | Piperacillin/tazobactam + Tobramycin |
|
||||
|:------------------------|:-------------------------------------|:-------------------------------------|
|
||||
| 70.2% (65-75.4%) | 93.6% (92-95.1%) | 89.9% (87-92.5%) |
|
||||
| 72.5% (67.7-77.1%) | 96% (95.1-96.9%) | 92.1% (89.3-94.4%) |
|
||||
|
||||
WISCA supports stratification by any clinical variable, so you can
|
||||
generate syndrome-specific or ward-specific coverage estimates:
|
||||
@@ -234,13 +234,19 @@ wisca(example_isolates,
|
||||
antimicrobials = c("TZP", "TZP+TOB", "TZP+GEN"),
|
||||
syndromic_group = "ward",
|
||||
minimum = 10) # Recommended threshold: >=30
|
||||
#> Warning: In 6 of 9 WISCA regimen(s) or group(s), the prior weight of 1 per pathogen accounts for more than
|
||||
#> 5% of the pathogen incidence (up to 13.3%), since many pathogens occur only rarely. Because WISCA
|
||||
#> gives every pathogen a prior weight of 1 (Bielicki et al., 2016), such rare pathogens can strongly
|
||||
#> lower the estimated coverage without reflecting the clinical expectation. Consider limiting the
|
||||
#> analysis to the most common pathogens with `top_n_microorganisms()`, or grouping rare pathogens
|
||||
#> before running WISCA.
|
||||
```
|
||||
|
||||
| Syndromic Group | Piperacillin/tazobactam | Piperacillin/tazobactam + Gentamicin | Piperacillin/tazobactam + Tobramycin |
|
||||
|:----------------|:------------------------|:-------------------------------------|:-------------------------------------|
|
||||
| Clinical | 74.8% (69-80.3%) | 93.7% (92.1-95.1%) | 90.4% (86.8-93.3%) |
|
||||
| ICU | 56.7% (48.2-65.4%) | 86.6% (83.3-89.6%) | 83.1% (78.2-87.5%) |
|
||||
| Outpatient | 57.6% (46-69%) | 76.4% (69.8-81.7%) | 68% (57.7-77.5%) |
|
||||
| Clinical | 78% (72.7-83.4%) | 96.3% (95.1-97.3%) | 93% (89.4-95.7%) |
|
||||
| ICU | 59.7% (53.5-67.9%) | 92.2% (90-94.3%) | 88.2% (83.4-92.2%) |
|
||||
| Outpatient | 67.9% (51.9-84.7%) | 92.8% (88.2-96.6%) | 79.8% (64.9-92.9%) |
|
||||
|
||||
**For AMR surveillance**, traditional antibiograms remain the right tool
|
||||
for tracking resistance per species over time:
|
||||
@@ -427,7 +433,7 @@ data analysis. This package can be used for:
|
||||
organisms (MDRO) ([tutorial](./reference/mdro.html))
|
||||
- Calculating (empirical) susceptibility of both mono therapy and
|
||||
combination therapies ([tutorial](./articles/AMR.html))
|
||||
- Apply AMR functions in predictive modelling
|
||||
- Applying AMR functions in predictive modelling
|
||||
([tutorial](./articles/AMR_with_tidymodels.html))
|
||||
- Getting properties for any microorganism (like Gram stain, species,
|
||||
genus or family) ([manual](./reference/mo_property.html))
|
||||
@@ -435,7 +441,8 @@ data analysis. This package can be used for:
|
||||
EARS-Net/ATC/LOINC/PubChem, defined daily dose or trade name)
|
||||
([manual](./reference/ab_property.html))
|
||||
- Plotting antimicrobial resistance ([tutorial](./articles/AMR.html))
|
||||
- Applying EUCAST expert rules ([manual](./reference/eucast_rules.html))
|
||||
- Applying EUCAST and CLSI interpretive rules
|
||||
([manual](./reference/interpretive_rules.html))
|
||||
- Getting SNOMED codes of a microorganism, or getting properties of a
|
||||
microorganism based on a SNOMED code
|
||||
([manual](./reference/mo_property.html))
|
||||
|
||||
@@ -268,6 +268,10 @@ It does so by combining two quantities: the relative incidence of each pathogen
|
||||
|
||||
\strong{Prior distributions:} Pathogen incidence uses a non-informative \eqn{Dirichlet(1, 1, \ldots, 1)} prior. Susceptibility proportions use the Jeffreys prior, \eqn{\beta(0.5, 0.5)}, except for pathogen-drug combinations with known intrinsic resistance, which use a strongly informative \eqn{\beta(1, 9999)} prior that forces near-zero susceptibility regardless of observed data. Intrinsic resistance is determined using the \link{intrinsic_resistant} data set, which is based on \href{https://www.eucast.org/bacteria/important-additional-information/expert-rules/}{'EUCAST Expected Resistant Phenotypes' v1.2} (2023).
|
||||
|
||||
\strong{Pathogens included in the model:} For each regimen (and group), only pathogens that occur in the data and have at least one test result for the regimen are included. Pathogens without any test result for a regimen do not contribute to its coverage, nor to its \code{n_total} in \code{long_numeric}. Use \code{\link[=retrieve_wisca_parameters]{retrieve_wisca_parameters()}} to get the parameters that were used.
|
||||
|
||||
\strong{Rare pathogens:} Since every pathogen receives a prior weight of 1 in the Dirichlet distribution, pathogens found only once or a few times are given relatively much weight, which pulls the incidence towards a uniform distribution. Their very uncertain susceptibility can then strongly lower the estimated coverage, without reflecting the clinical expectation. A warning is therefore given when this prior weight accounts for more than 5\% of the pathogen incidence in any regimen (or group). Consider limiting the analysis to the most common pathogens using \code{\link[=top_n_microorganisms]{top_n_microorganisms()}}, or grouping rare pathogens before running WISCA.
|
||||
|
||||
\strong{Interpreting the output:} Overlapping credible intervals between regimens indicate no significant difference in coverage; if a narrower-spectrum regimen overlaps with a broader one, the narrower-spectrum option may be preferred on stewardship grounds. Non-overlapping intervals indicate a clinically meaningful difference. For small sample sizes, consider pooling data from multiple sites to improve precision, provided pathogen distributions are sufficiently similar (Bielicki \emph{et al.}, 2016).
|
||||
|
||||
For the full mathematical derivation and worked examples, see the \href{https://amr-for-r.org/articles/WISCA.html}{WISCA vignette}.
|
||||
|
||||
@@ -47,6 +47,35 @@ test_that("test-antibiogram.R", {
|
||||
wisca_known <- suppressMessages(wisca(example_isolates[example_isolates$mo != "UNKNOWN", ], antimicrobials = "GEN", simulations = 100, info = FALSE, guideline = "EUCAST"))
|
||||
expect_true("UNKNOWN" %in% attributes(wisca_all)$wisca_parameters$mo)
|
||||
expect_false(identical(as.data.frame(wisca_all), as.data.frame(wisca_known)))
|
||||
|
||||
# WISCA only includes pathogens that occur in a group and were tested for the regimen (#337)
|
||||
set.seed(337)
|
||||
wisca_data <- example_isolates[mo_domain(example_isolates$mo) == "Bacteria", ]
|
||||
wisca_data$grp <- sample(c("A", "B", "C"), nrow(wisca_data), replace = TRUE)
|
||||
wisca_strat <- suppressWarnings(wisca(wisca_data, antimicrobials = "MEM", syndromic_group = "grp", simulations = 5000, info = FALSE))
|
||||
wisca_params <- retrieve_wisca_parameters(wisca_strat)
|
||||
# no absent or untested pathogens
|
||||
expect_true(all(wisca_params$n_total > 0 & wisca_params$n_tested > 0))
|
||||
# each stratum equals the unstratified WISCA of its subset
|
||||
for (g in c("A", "B", "C")) {
|
||||
wisca_subset <- suppressWarnings(wisca(wisca_data[wisca_data$grp == g, ], antimicrobials = "MEM", simulations = 5000, info = FALSE))
|
||||
cov_strat <- attributes(wisca_strat)$long_numeric$coverage[attributes(wisca_strat)$long_numeric$syndromic_group == g]
|
||||
expect_true(abs(cov_strat - attributes(wisca_subset)$long_numeric$coverage) < 0.01)
|
||||
expect_identical(attributes(wisca_strat)$long_numeric$n_total[attributes(wisca_strat)$long_numeric$syndromic_group == g], attributes(wisca_subset)$long_numeric$n_total)
|
||||
}
|
||||
# untested pathogens have no influence
|
||||
wisca_untested <- wisca_data[1:20, ]
|
||||
wisca_untested$mo <- as.mo("Brucella melitensis")
|
||||
wisca_untested$MEM <- NA_sir_
|
||||
set.seed(337)
|
||||
wisca_without <- suppressWarnings(wisca(wisca_data, antimicrobials = "MEM", simulations = 1000, info = FALSE))
|
||||
set.seed(337)
|
||||
wisca_with <- suppressWarnings(wisca(rbind(wisca_data, wisca_untested), antimicrobials = "MEM", simulations = 1000, info = FALSE))
|
||||
expect_identical(attributes(wisca_with)$long_numeric, attributes(wisca_without)$long_numeric)
|
||||
expect_false("B_BRCLL_MLTN" %in% retrieve_wisca_parameters(wisca_with)$mo)
|
||||
# warning for a large prior weight of rare pathogens
|
||||
expect_warning(wisca(wisca_data, antimicrobials = "MEM", syndromic_group = "grp", simulations = 100, info = FALSE), "prior weight")
|
||||
expect_warning(wisca(example_isolates, antimicrobials = "MEM", simulations = 100, info = FALSE), NA)
|
||||
if (AMR:::pkg_is_available("dplyr", min_version = "1.0.0", also_load = TRUE)) {
|
||||
ab_grouped <- tryCatch(
|
||||
# pathogen names are translated, so test in English, which must also hold on systems in another language
|
||||
@@ -230,7 +259,19 @@ test_that("test-antibiogram.R", {
|
||||
expect_equal(nrow(ab_grp_seq), nrow(ab_grp_par))
|
||||
}
|
||||
|
||||
# 6. parallel = TRUE without a plan raises an informative error
|
||||
# 6. parallel and sequential use the same WISCA parameters and give the same coverage (#337)
|
||||
set.seed(337)
|
||||
wisca_data <- example_isolates[mo_domain(example_isolates$mo) == "Bacteria", ]
|
||||
wisca_data$grp <- sample(c("A", "B", "C"), nrow(wisca_data), replace = TRUE)
|
||||
future::plan(future::sequential)
|
||||
wisca_strat_seq <- suppressWarnings(wisca(wisca_data, antimicrobials = "MEM", syndromic_group = "grp", simulations = 5000, info = FALSE))
|
||||
future::plan(future::multicore)
|
||||
wisca_strat_par <- suppressWarnings(wisca(wisca_data, antimicrobials = "MEM", syndromic_group = "grp", simulations = 5000, parallel = TRUE, info = FALSE))
|
||||
expect_identical(retrieve_wisca_parameters(wisca_strat_par), retrieve_wisca_parameters(wisca_strat_seq))
|
||||
expect_true(all(retrieve_wisca_parameters(wisca_strat_par)$n_tested > 0))
|
||||
expect_true(all(abs(attributes(wisca_strat_par)$long_numeric$coverage - attributes(wisca_strat_seq)$long_numeric$coverage) < 0.01))
|
||||
|
||||
# 7. parallel = TRUE without a plan raises an informative error
|
||||
future::plan(future::sequential)
|
||||
expect_error(
|
||||
suppressWarnings(wisca(example_isolates, antimicrobials = "TZP", parallel = TRUE, info = FALSE)),
|
||||
|
||||
Reference in new issue
Block a user