Affiliations
Sydney School of Veterinary Science, Faculty of Science, The
University of Sydney, Camden, New South Wales 2570, Australia
library(tidyverse)
library(purrr)
library(readxl)
library(janitor)
#library(table1)
# install.packages("furniture")
Import enrollment data
enrollment <- read_excel("Bulk milk AMR data entry.xlsx",
col_types = c("numeric", "text", "numeric",
"numeric", "numeric", "numeric",
"numeric", "text", "date", "date",
"date", "date", "date", "date", "text",
"text", "text", "text", "text", "text",
"numeric", "text", "text", "numeric",
"numeric", "numeric", "numeric",
"text", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric","text","text","text","text"))
enrollment <- enrollment %>% clean_names() %>% select(ref,
sample_1_date_of_collection,
sample_1_date_of_freezing,
sample_1_date_of_culture,
sample_2_date_of_collection,
sample_2_date_of_freezing,
sample_2_date_of_culture,
region,
breed,
ave_milk_production,
most_recent_bulk_milk_scc,
pasture_y_n,
mixed_ratio_y_n,
no_lactating_cows,
outlet_1_or_sample_port_2
)
enrollment$farm_id <- enrollment$ref
enrollment$sample_1_freeze_lag <- difftime(enrollment$sample_1_date_of_freezing,enrollment$sample_1_date_of_collection, units = "days")
enrollment$sample_1_culture_lag <- difftime(enrollment$sample_1_date_of_culture,enrollment$sample_1_date_of_freezing, units = "days")
enrollment$sample_2_freeze_lag <- difftime(enrollment$sample_2_date_of_freezing,enrollment$sample_2_date_of_collection, units = "days")
enrollment$sample_2_culture_lag <- difftime(enrollment$sample_2_date_of_culture,enrollment$sample_2_date_of_freezing, units = "days")
print(summarytools::dfSummary(enrollment,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")
table1::table1(~ breed + ave_milk_production + most_recent_bulk_milk_scc + factor(pasture_y_n) + factor(mixed_ratio_y_n) + no_lactating_cows + factor(outlet_1_or_sample_port_2), data=enrollment)
|
Overall
(N=40) |
| breed |
|
| Ayrshire |
1 (2.5%) |
| Brown Swiss |
2 (5.0%) |
| Guernsey |
1 (2.5%) |
| Holstein Friesian |
30 (75.0%) |
| Holstein Friesian & Jersey |
6 (15.0%) |
| ave_milk_production |
|
| Mean (SD) |
24.0 (5.54) |
| Median [Min, Max] |
24.0 [15.0, 40.0] |
| most_recent_bulk_milk_scc |
|
| Mean (SD) |
208 (67.5) |
| Median [Min, Max] |
200 [90.0, 426] |
| factor(pasture_y_n) |
|
| 0 |
5 (12.5%) |
| 1 |
35 (87.5%) |
| factor(mixed_ratio_y_n) |
|
| 0 |
23 (57.5%) |
| 1 |
17 (42.5%) |
| no_lactating_cows |
|
| Mean (SD) |
520 (1220) |
| Median [Min, Max] |
270 [87.0, 7800] |
| factor(outlet_1_or_sample_port_2) |
|
| 0 |
34 (85.0%) |
| 1 |
3 (7.5%) |
| 2 |
3 (7.5%) |
table1::table1(~ as.numeric(sample_1_freeze_lag) + as.numeric(sample_1_culture_lag) + as.numeric(sample_2_freeze_lag) + as.numeric(sample_2_culture_lag) , data=enrollment)
|
Overall
(N=40) |
| as.numeric(sample_1_freeze_lag) |
|
| Mean (SD) |
1.30 (1.40) |
| Median [Min, Max] |
1.00 [0, 4.00] |
| as.numeric(sample_1_culture_lag) |
|
| Mean (SD) |
54.8 (14.2) |
| Median [Min, Max] |
56.0 [16.0, 75.0] |
| as.numeric(sample_2_freeze_lag) |
|
| Mean (SD) |
4.03 (16.1) |
| Median [Min, Max] |
1.00 [0, 103] |
| as.numeric(sample_2_culture_lag) |
|
| Mean (SD) |
60.5 (32.7) |
| Median [Min, Max] |
53.5 [11.0, 200] |
result <- read_excel("Bulk milk AMR data entry.xlsx",
sheet = "results")
result <- result %>% clean_names %>% select(farm_id,round,
esbl,mac,mac_lact,esbl_1,esbl_id_1,
staph_aureus,staph24_1,staph24_id_1,staph24_2,staph24_id_2,
mrsa,mrsa_1,mrsa_1_id,mrsa_2,mrsa_2_id,
strep_ag,edwards,strep_beta,strep_beta_id,
vre,vre_1,vre_1_id,vre_2,vre_2_id
)
staph_aureus <- result %>% group_by(farm_id) %>% summarise(
staph_aureus_any = sum(staph_aureus==1)
)
#
# staph_aureus <- merge(staph_aureus,enrollment %>% select(farm_id,most_recent_bulk_milk_scc))
# staph_aureus$staph_aureus <- 0
# staph_aureus$staph_aureus[staph_aureus$staph_aureus_any >= 1] <- 1
#
# staph_aureus %>% group_by(staph_aureus) %>% summarise(
# ave_scc <- mean(most_recent_bulk_milk_scc, na.rm=T)
# )
print(summarytools::dfSummary(result,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")
print(summarytools::dfSummary(staph_aureus,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")
Data Frame Summary
staph_aureus
Dimensions: 40 x 2
Duplicates: 0
| No |
Variable |
Stats / Values |
Freqs (% of Valid) |
Graph |
Missing |
| 1 |
farm_id
[numeric] |
| Mean (sd) : 20.5 (11.7) | | min ≤ med ≤ max: | | 1 ≤ 20.5 ≤ 40 | | IQR (CV) : 19.5 (0.6) |
|
40 distinct values |
 |
0
(0.0%) |
| 2 |
staph_aureus_any
[integer] |
| Mean (sd) : 0.4 (0.6) | | min ≤ med ≤ max: | | 0 ≤ 0 ≤ 2 | | IQR (CV) : 1 (1.7) |
|
| 0 | : | 28 | ( | 70.0% | ) | | 1 | : | 9 | ( | 22.5% | ) | | 2 | : | 3 | ( | 7.5% | ) |
|
 |
0
(0.0%) |
Generated by summarytools 1.0.1 (R version 4.2.1)
2023-05-01
table1::table1(~ esbl_id_1 + staph24_id_1 + staph24_id_2 + mrsa_1_id + mrsa_2_id + vre_1_id | round,data=result)
|
1
(N=40) |
2
(N=40) |
Overall
(N=80) |
| esbl_id_1 |
|
|
|
| No ID |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Pseudomonas aeruginosa |
5 (12.5%) |
2 (5.0%) |
7 (8.8%) |
| Pseudomonas protegens |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Missing |
33 (82.5%) |
38 (95.0%) |
71 (88.8%) |
| staph24_id_1 |
|
|
|
| Catalase negative |
11 (27.5%) |
13 (32.5%) |
24 (30.0%) |
| No ID |
3 (7.5%) |
2 (5.0%) |
5 (6.3%) |
| Staphylococcus aureus |
5 (12.5%) |
9 (22.5%) |
14 (17.5%) |
| Staphylococcus chromogenes |
2 (5.0%) |
3 (7.5%) |
5 (6.3%) |
| Staphylococcus haemolyticus |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Staphylococcus simulans |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Stenotrophomonas maltophilia |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Streptococcus uberis |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Corynebacterium provencense |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Diutina rugosa |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Enterococcus saccharolyticus |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Kurthia gibsonii |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Missing |
15 (37.5%) |
9 (22.5%) |
24 (30.0%) |
| staph24_id_2 |
|
|
|
| Catalase negative |
2 (5.0%) |
0 (0%) |
2 (2.5%) |
| Staphylococcus aureus |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Missing |
37 (92.5%) |
40 (100%) |
77 (96.3%) |
| mrsa_1_id |
|
|
|
| Arthrobacter woluwensis |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Catalase negative |
10 (25.0%) |
4 (10.0%) |
14 (17.5%) |
| Cellulomonas sp. |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| No ID |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Staphylococcus pasteuri |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Lactococcus raffinolactis |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Macrococcus caseolyticus |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Missing |
26 (65.0%) |
34 (85.0%) |
60 (75.0%) |
| mrsa_2_id |
|
|
|
| Catalase negative |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Microbacterium sp. |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| No ID |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Missing |
37 (92.5%) |
40 (100%) |
77 (96.3%) |
| vre_1_id |
|
|
|
| Catalase positive |
6 (15.0%) |
4 (10.0%) |
10 (12.5%) |
| Sphingobacterium multivorum |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Sphingobacterium sp. |
1 (2.5%) |
0 (0%) |
1 (1.3%) |
| Sphingomonas sp. |
0 (0%) |
1 (2.5%) |
1 (1.3%) |
| Missing |
32 (80.0%) |
35 (87.5%) |
67 (83.8%) |
result %>% select(staph24_id_1,staph24_id_2) %>% pivot_longer(cols = starts_with("staph24"),names_to = "isolate",values_to = "id") %>% filter(!is.na(id)) %>% tabyl(id) %>% arrange(desc(percent))
## id n percent
## Catalase negative 26 0.44067797
## Staphylococcus aureus 15 0.25423729
## No ID 5 0.08474576
## Staphylococcus chromogenes 5 0.08474576
## Corynebacterium provencense 1 0.01694915
## Diutina rugosa 1 0.01694915
## Enterococcus saccharolyticus 1 0.01694915
## Kurthia gibsonii 1 0.01694915
## Staphylococcus haemolyticus 1 0.01694915
## Staphylococcus simulans 1 0.01694915
## Stenotrophomonas maltophilia 1 0.01694915
## Streptococcus uberis 1 0.01694915
result %>% select(mrsa_1_id,mrsa_2_id) %>% pivot_longer(cols = starts_with("mrsa"),names_to = "isolate",values_to = "id") %>% filter(!is.na(id)) %>% tabyl(id) %>% arrange(desc(percent))
## id n percent
## Catalase negative 15 0.65217391
## No ID 2 0.08695652
## Arthrobacter woluwensis 1 0.04347826
## Cellulomonas sp. 1 0.04347826
## Lactococcus raffinolactis 1 0.04347826
## Macrococcus caseolyticus 1 0.04347826
## Microbacterium sp. 1 0.04347826
## Staphylococcus pasteuri 1 0.04347826
result %>% select(vre_1_id,vre_2_id) %>% pivot_longer(cols = starts_with("vre"),names_to = "isolate",values_to = "id") %>% filter(!is.na(id)) %>% tabyl(id) %>% arrange(desc(percent))
## id n percent
## Catalase positive 11 0.78571429
## Sphingobacterium multivorum 1 0.07142857
## Sphingobacterium sp. 1 0.07142857
## Sphingomonas sp. 1 0.07142857
# library(furniture)
# table1 <- table1(enrollment,second = c("no_lactating_cows","most_recent_bulk_milk_scc","ave_milk_production"),
# region,breed,no_lactating_cows,ave_milk_production,most_recent_bulk_milk_scc,factor(pasture_y_n))
#
# #write.csv(table1,"table1.csv")
# table1
Sensitivity analysis
The objective of our sensitivity analysis was to estimate the
probability of all sampled herds testing negative under different herd
prevalence and test sensitivity scenarios. To do this, we simulated
populations of NSW dairy herds (n = 494). For each iteration (n =
40,000) a value was chosen for herd prevalence (either 5%, 10%, 20%, or
30%) and test sensitivity (randomly selected from a uniform distribution
between 0 and 1). Simulated herds were classified as being positive or
negative for the target organism (e.g., MRSA), to match the chosen
overall herd prevalence. Forty herds were then randomly selected from
each simulated population and 2 Bernoulli trials (probability = test
sensitivity) were used to simulate testing of each herd. It was assumed
that the test had 100% specificity because species identity and AMR
status would have been confirmed with additional laboratory methods. If
a simulated herd had at least 1 test positive result, it was determined
to have been detected by the testing regime. A regression line
(logistic) was used to plot the relationship between test sensitivity
and probability of all herds testing negative.
study_sim <- function(herd_prev_pop) {
sens <- runif(1, min = 0, max = 1)
n_infected <- round(herd_prev_pop * 494, 0)
pop <- tibble(
herd = 1:494
)
pop <- pop %>% mutate(
mrsa_true = case_when(
herd <= n_infected ~ 1,
herd > n_infected ~ 0
),
mrsa_apparent_1 = case_when(
mrsa_true == 0 ~ 0,
mrsa_true == 1 ~ rbinom(494, 1, sens)
),
mrsa_apparent_2 = case_when(
mrsa_true == 0 ~ 0,
mrsa_true == 1 ~ rbinom(494, 1, sens)
)
)
sample <- sample_n(pop, 40)
output <- tibble(
sensitivity = sens,
herd_prev_true = herd_prev_pop,
test_1 = sample$mrsa_apparent_1,
test_2 = sample$mrsa_apparent_2
)
output %>% summarise(
sensitivity = sens,
herd_prev_true = herd_prev_pop,
test_pos = sum(test_1 == 1 | test_2 == 1)
)
}
set.seed(1234)
values <- c(0.05, 0.1, 0.2, 0.3)
result_list <- map_dfr(values, function(x) {
map_dfr(rep(1:10000, simplify = FALSE), ~ study_sim(herd_prev_pop = x), .id = "replicate")
}, .id = "value")
total <- result_list %>%
mutate(
sensitivity = as.numeric(sensitivity),
test_pos = as.numeric(test_pos),
herd_prev_true = as.numeric(herd_prev_true) %>% as.factor(),
at_least_one_herd = ifelse(test_pos > 0, 1, 0),
all_herds_negative = ifelse(at_least_one_herd == 1, 0, 1)
)
plot <- ggplot(total) + coord_cartesian(ylim = (c(0,1))) +
aes(x=sensitivity, y=all_herds_negative, group=herd_prev_true, colour=herd_prev_true) +
labs(y = "Probability of all herds testing negative", x = "Test Sensitivity") +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE) +
scale_y_continuous(breaks=seq(0, 1, .1)) +
scale_x_continuous(breaks=seq(0, 1, .1)) +
theme(panel.border = element_rect(colour = "black",fill=NA, size=1),
axis.text=element_text(size=12,family="Times"),
axis.title=element_text(size=14,face="bold",family="Times"),
panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),
legend.position=c(0.9,0.8),
legend.text = element_text(colour="black", size=12,face="bold",family="Times")) +
scale_color_discrete(name = "Prevalence", labels = c("5%", "10%", "20%","30%"))
plot
