We would like to thank the participating farmers and laboratory staff for their essential roles in the provision and testing of samples.
This study was conducted as part of Food Agility project FA 078 – On-farm mastitis, jointly funded by Food Agility CRC, Dairy Australia, Coles Sustainable Dairy Development Group, Charles Sturt University, University of Technology Sydney, and the University of Sydney.
The primary objective of this study was to evaluate the diagnostic performance of 5 point of care tests available in Australia for guiding mastitis therapy. A secondary objective was to describe the pathogen profiles of mastitis-causing organisms in cows managed in barns (“intensive”) and on pasture (“non-intensive”).
library(tidyverse)
library(janitor)
library(epiR)
library(DT)
library(ggcharts)
all_samples <- readRDS("~/Dropbox/R backup/2022/da_mastitis_tests/all_samples.rds")
data <- all_samples
data$farm_id <- NULL
data$farm_code_other <- NULL
data$sample_date <- as.Date(data$sample_date)
data$test_date <- as.Date(data$test_date)
data$test_lag <- as.numeric(data$test_date - data$sample_date)
#print(summarytools::dfSummary(data,valid.col=FALSE, graph.magnif=0.8, style="grid"),method = "render")A total of 641 samples from 30 farms were tested during the study period. The number of samples tested from each farm ranged from 1 to 82 (median = 15). Samples were submitted from intensive (n = 207 samples from 6 farms) and non-intensive herds (n = 434 samples from 24 farms) in New South Wales (n = 360 samples from 15 farms) and Victoria (n = 281 samples from 15 farms). Contamination was identified in 149 samples (23.2%). For farms that submitted 10 or more samples, the herd-level contamination prevalence ranged from 0.0 to 47.4% (median = 22.2%; Figure 1). The percentage of samples testing positive for each of the diagnoses of interest are shown in Figure 2. After excluding contaminated samples, the prevalence of diagnoses of interest for samples from intensive herds in descending order were no growth (31.7%), Klebsiella spp. (28.1%), E. coli (15.0%), Strep. uberis (8.4%), Strep. dysgalactiae (2.4%), and Staph. aureus (0.0%). Prevalences in non-contaminated samples from non-intensive herds were Strep. uberis (35.0%), no growth (26.9%), E. coli (13.3%), Strep. dysgalactiae (5.0%), Staph. aureus (4.3%), and Klebsiella spp. (1.5%).
data %>% tabyl(contamination)| contamination | n | percent |
|---|---|---|
| 0 | 492 | 0.7675507 |
| 1 | 149 | 0.2324493 |
data %>% tabyl(contamination,system)| contamination | Intensive | Non-intensive | NA_ |
|---|---|---|---|
| 0 | 167 | 323 | 2 |
| 1 | 40 | 106 | 3 |
#a <- data$test_date %>% as.data.frame()
#b <- data$sample_date %>% as.data.frame()herd_contamination <- data %>% tabyl(farm_number,contamination) %>% as_tibble
herd_contamination <- herd_contamination %>% mutate(
samples = `0` + `1`,
contaminated = `1` / samples,
sample_gt_9 = case_when(
samples > 9 ~ 1,
T ~ 0
)
) %>% select(samples,contaminated,sample_gt_9)
herd_contamination %>% mutate(contaminated = round(contaminated,3)) %>% datatable(filter = "top")contaminated <- herd_contamination %>% filter(sample_gt_9 == 1) %>% select(contaminated)
contaminated$contaminated %>% median[1] 0.2222222ggplot(herd_contamination %>% filter(sample_gt_9 == 1), aes(x = contaminated)) +
geom_histogram(fill = "steelblue", color = "white",binwidth = 0.1) +
scale_y_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
scale_x_continuous(breaks = c(0, 0.1, 0.2, 0.3, 0.4, 0.5)) +
theme(panel.background = element_blank(),
panel.border = element_rect(color = "black", fill = NA, size = 1)) +
labs(x = "Proportion of samples that were contaminated (3 of more morphotypes)", y = "Frequency of herds") +
ggtitle("Contamination in herds with 10 or more samples")
ggsave("contamination.tiff", width = 6, height = 6, dpi = 200)This table shows the frequency of culture results from the 492 samples that were not contaminated.
freq <- rbind(data$iso1_dx,
data$iso2_dx) %>% as.vector %>% tabyl %>% tibble
valid <- sum((data$contamination==0) == T)
freq$pathogen <- freq$.
freq <- freq %>% mutate(
percent = round(n / valid,3)*100
) %>% filter(!is.na(pathogen))
freq <- freq %>% select(pathogen,n,percent)
freq %>% datatable(filter = "top")# Pathogen profiles
result_out <- data %>% filter(contamination==0) %>% group_by(system) %>% summarise(
`No growth` = sum(no_growth==1),
`E. coli` = sum(e_coli==1),
`Klebsiella spp.` = sum(klebsiella==1),
`Staph. aureus` = sum(staph_aureus==1),
`Strep. dysgalactiae` = sum(strep_dys==1),
`Strep. uberis` = sum(strep_uberis==1)) %>% pivot_longer(`No growth`:`Strep. uberis`,
names_to = "Result",
values_to = "count")
n_intensive_non_contaminated <- data %>% filter(system=="Intensive" & contamination==0) %>% count %>% as.numeric
n_non_intensive_non_contaminated <- data %>% filter(system=="Non-intensive" & contamination==0) %>% count %>% as.numeric
result_out <- result_out %>% mutate(
Percent = case_when(
system=="Intensive" ~ count / n_intensive_non_contaminated,
system=="Non-intensive" ~ count / n_non_intensive_non_contaminated
)
)
#intensive
profile_intensive <- ggplot(result_out %>% filter(system=="Intensive")) +
aes(y = Percent, reorder(Result,-Percent)) +
geom_bar(position = "dodge", stat="identity",color='skyblue',fill='steelblue') +
scale_y_continuous(breaks=seq(0, 0.5, 0.1),labels = scales::percent, limits=c(0, 0.4)) +
ggtitle("Intensive") +
labs(y = "Percentage of cows", x = "") +
scale_fill_brewer(palette="Set1") +
# geom_text(aes(x = Test, y = Percent), size=3.5, nudge_y = .02 ,colour='black') +
labs(y = "", x = "") +
#geom_text(aes(x = culture_result, y = Percent, fill = Test, label=Percent), size=3.5, nudge_y = .02 ,colour='black')
theme(panel.border = element_rect(colour = "black", fill=NA,size=1),
title=element_text(size=12,face="bold",family="Times"),
axis.text=element_text(size=12,family="Times"),
axis.title=element_text(size=11,face="bold",family="Times"),
axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.2, linetype = 'solid', colour = "grey"),
panel.grid.minor = element_line(size = 0.2, linetype = 'solid',colour = "grey"),
legend.text = element_text(colour="black",size=8,face="bold",family="Times"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
profile_non_intensive <- ggplot(result_out %>% filter(system=="Non-intensive")) +
aes(y = Percent, x = reorder(Result,-Percent)) +
geom_bar(position = "dodge", stat="identity",color='skyblue',fill='steelblue') +
scale_y_continuous(breaks=seq(0, 1, 0.1),labels = scales::percent, limits=c(0, 0.4)) +
ggtitle("Non-intensive") +
labs(y = "Percentage of cows", x = "") +
scale_fill_brewer(palette="Set1") +
# geom_text(aes(x = Test, y = Percent), size=3.5, nudge_y = .02 ,colour='black') +
labs(y = "", x = "") +
#geom_text(aes(x = culture_result, y = Percent, fill = Test, label=Percent), size=3.5, nudge_y = .02 ,colour='black')
theme(panel.border = element_rect(colour = "black", fill=NA,size=1),
title=element_text(size=12,face="bold",family="Times"),
axis.text=element_text(size=12,family="Times"),
axis.title=element_text(size=11,face="bold",family="Times"),
axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.2, linetype = 'solid', colour = "grey"),
panel.grid.minor = element_line(size = 0.2, linetype = 'solid',colour = "grey"),
legend.text = element_text(colour="black",size=8,face="bold",family="Times"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())This table shows the pathogen profiles cows in intensive and non-intensive farming systems.
#result_out
# ggplot(data = result_out %>% filter(!is.na(system))) +
# facet_wrap(~system,ncol=2) +
# aes(y=Percent,x = Result) +
# geom_bar(position = "dodge", stat="identity",color='skyblue',fill='steelblue') +
# theme(legend.position="bottom",
# legend.direction="horizontal",
# panel.background = element_blank(),
# legend.title = element_blank(),
# panel.border = element_rect(color = "grey", fill = NA),
# text = element_text(family = "Times New Roman"),
# axis.ticks.x = element_blank(),
# panel.grid.major.y = element_line(color = "gray80", linetype = "dashed", size = 0.5),
# panel.grid.minor.y = element_line(color = "gray90", linetype = "dotted", size = 0.5)) +
# scale_y_continuous(limits = c(0, .5), breaks = seq(0, 0.5,0.1), minor_breaks = seq(0, 0.5, 0.1)) +
# labs(x="", y="Percentage") +
# ggtitle("Non-intensive herds")ggplot(data = result_out %>% filter(!is.na(system))) +
aes(y = Percent, x = reorder(Result, -Percent), fill = system) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = paste0(round(Percent * 100, 1), "%")), position = position_dodge(width = 0.9), vjust = -0.5, size = 3) +
theme(
legend.position = c(0.8, 0.7),
legend.direction = "horizontal",
panel.background = element_blank(),
legend.title = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman"), # Set font to Times New Roman
axis.ticks.x = element_blank()
) +
scale_y_continuous(limits = c(0, 0.4),
breaks = seq(0, 0.4, 0.1),
minor_breaks = seq(0, 0.4, 0.1),
labels = scales::percent) +
labs(x = "", y = "Percentage of cases") +
ggtitle("Pathogen profile")
ggsave("pathogen_profile.tiff", width = 6, height = 4, dpi = 200)This table shows the sample data as the bar chart above.
result_out %>% filter(!is.na(system)) %>% mutate(
Percent = round(Percent* 100,1)
) %>% datatable(filter = "top")These tables are the same data again, but in descending order.
profile_non_intensive
profile_intensive
ab_use_intensive<- data %>% filter(contamination==0 & system=="Intensive") %>% summarise(
culture_result = "Non-contaminated samples",
sample_n= sum(contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
ab_use_non_intensive<- data %>% filter(contamination==0 & system=="Non-intensive") %>% summarise(
culture_result = "Non-contaminated samples",
sample_n= sum(contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
# Antibiotic use stratified by system
ab_use_intensive$System <- "Intensive"
ab_use_non_intensive$System <- "Non-intensive"
ab_use <- rbind(ab_use_intensive,
ab_use_non_intensive) %>% select(System,Biplate:Mastatest) %>% pivot_longer(Biplate:Mastatest,
values_to = "Percent",
names_to = "Test")
ab_use$Percent <- ab_use$Percent / 100
plot2 <- ggplot(ab_use) +
aes(y = Percent, x = System) +
geom_bar(aes(fill = Test), position = "dodge", stat="identity") +
scale_y_continuous(breaks=seq(0, 1, 0.2),labels = scales::percent) +
ggtitle(paste0("Probability of antibiotic treatment")) +
labs(y = "Percentage of cows", x = "") +
scale_fill_brewer(palette="Set1") +
# geom_text(aes(x = Test, y = Percent), size=3.5, nudge_y = .02 ,colour='black') +
labs(y = "", x = "") +
#geom_text(aes(x = culture_result, y = Percent, fill = Test, label=Percent), size=3.5, nudge_y = .02 ,colour='black')
theme(panel.border = element_rect(colour = "black", fill=NA,size=1),
title=element_text(size=12,face="bold",family="Times"),
axis.text=element_text(size=12,family="Times"),
axis.title=element_text(size=11,face="bold",family="Times"),
axis.text.x = element_text(size = 10),
panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.2, linetype = 'solid', colour = "grey"),
panel.grid.minor = element_line(size = 0.2, linetype = 'solid',colour = "grey"),
legend.text = element_text(colour="black",size=8,face="bold",family="Times"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())The figure below shows the percentage of cases identified as antibiotic treatment candidates using 5 point of care testing systems. One clear finding here is that a lower proportion cases in the intensive farming system are identified as antibiotic treatment candidates.
plot2
ab_use_intensive_stack <- data %>% filter(contamination==0 & system=="Intensive") %>% summarise(
culture_result = "Non-contaminated samples",
sample_n= sum(contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate_tx1_dz1 = sum(biplate_tx_gram_pos==1 & gram_pos==1,na.rm=T),
biplate_tx0_dz1 = sum(biplate_tx_gram_pos==0 & gram_pos==1,na.rm=T),
biplate_tx1_dz0 = sum(biplate_tx_gram_pos==1 & gram_pos==0,na.rm=T),
biplate_tx0_dz0 = sum(biplate_tx_gram_pos==0 & gram_pos==0,na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast_tx1_dz1 = sum(accumast_tx_gram_pos==1 & gram_pos==1,na.rm=T),
accumast_tx0_dz1 = sum(accumast_tx_gram_pos==0 & gram_pos==1,na.rm=T),
accumast_tx1_dz0 = sum(accumast_tx_gram_pos==1 & gram_pos==0,na.rm=T),
accumast_tx0_dz0 = sum(accumast_tx_gram_pos==0 & gram_pos==0,na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup_tx1_dz1 = sum(checkup_tx_gram_pos==1 & gram_pos==1,na.rm=T),
checkup_tx0_dz1 = sum(checkup_tx_gram_pos==0 & gram_pos==1,na.rm=T),
checkup_tx1_dz0 = sum(checkup_tx_gram_pos==1 & gram_pos==0,na.rm=T),
checkup_tx0_dz0 = sum(checkup_tx_gram_pos==0 & gram_pos==0,na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm_tx1_dz1 = sum(petrifilm_tx_gram_pos==1 & gram_pos==1,na.rm=T),
petrifilm_tx0_dz1 = sum(petrifilm_tx_gram_pos==0 & gram_pos==1,na.rm=T),
petrifilm_tx1_dz0 = sum(petrifilm_tx_gram_pos==1 & gram_pos==0,na.rm=T),
petrifilm_tx0_dz0 = sum(petrifilm_tx_gram_pos==0 & gram_pos==0,na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest_tx1_dz1 = sum(mastatest_tx_gram_pos==1 & gram_pos==1,na.rm=T),
mastatest_tx0_dz1 = sum(mastatest_tx_gram_pos==0 & gram_pos==1,na.rm=T),
mastatest_tx1_dz0 = sum(mastatest_tx_gram_pos==1 & gram_pos==0,na.rm=T),
mastatest_tx0_dz0 = sum(mastatest_tx_gram_pos==0 & gram_pos==0,na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
biplate_treat_true = round(100*biplate_tx1_dz1/(biplate_tx0+biplate_tx1),1),
biplate_treat_false = round(100*biplate_tx1_dz0/(biplate_tx0+biplate_tx1),1),
biplate_leave_true = round(100*biplate_tx0_dz0/(biplate_tx0+biplate_tx1),1),
biplate_leave_false = round(100*biplate_tx0_dz1/(biplate_tx0+biplate_tx1),1),
accumast_treat_true = round(100*accumast_tx1_dz1/(accumast_tx0+accumast_tx1),1),
accumast_treat_false = round(100*accumast_tx1_dz0/(accumast_tx0+accumast_tx1),1),
accumast_leave_true = round(100*accumast_tx0_dz0/(accumast_tx0+accumast_tx1),1),
accumast_leave_false = round(100*accumast_tx0_dz1/(accumast_tx0+accumast_tx1),1),
checkup_treat_true = round(100*checkup_tx1_dz1/(checkup_tx0+checkup_tx1),1),
checkup_treat_false = round(100*checkup_tx1_dz0/(checkup_tx0+checkup_tx1),1),
checkup_leave_true = round(100*checkup_tx0_dz0/(checkup_tx0+checkup_tx1),1),
checkup_leave_false = round(100*checkup_tx0_dz1/(checkup_tx0+checkup_tx1),1),
petrifilm_treat_true = round(100*petrifilm_tx1_dz1/(petrifilm_tx0+petrifilm_tx1),1),
petrifilm_treat_false = round(100*petrifilm_tx1_dz0/(petrifilm_tx0+petrifilm_tx1),1),
petrifilm_leave_true = round(100*petrifilm_tx0_dz0/(petrifilm_tx0+petrifilm_tx1),1),
petrifilm_leave_false = round(100*petrifilm_tx0_dz1/(petrifilm_tx0+petrifilm_tx1),1),
mastatest_treat_true = round(100*mastatest_tx1_dz1/(mastatest_tx0+mastatest_tx1),1),
mastatest_treat_false = round(100*mastatest_tx1_dz0/(mastatest_tx0+mastatest_tx1),1),
mastatest_leave_true = round(100*mastatest_tx0_dz0/(mastatest_tx0+mastatest_tx1),1),
mastatest_leave_false = round(100*mastatest_tx0_dz1/(mastatest_tx0+mastatest_tx1),1),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),
biplate_ppv = biplate_treat_true / (biplate_treat_true + biplate_treat_false),
biplate_npv = biplate_leave_true / (biplate_leave_true + biplate_leave_false),
accumast_ppv = accumast_treat_true / (accumast_treat_true + accumast_treat_false),
accumast_npv = accumast_leave_true / (accumast_leave_true + accumast_leave_false),
checkup_ppv = checkup_treat_true / (checkup_treat_true + checkup_treat_false),
checkup_npv = checkup_leave_true / (checkup_leave_true + checkup_leave_false),
petrifilm_ppv = petrifilm_treat_true / (petrifilm_treat_true + petrifilm_treat_false),
petrifilm_npv = petrifilm_leave_true / (petrifilm_leave_true + petrifilm_leave_false),
mastatest_ppv = mastatest_treat_true / (mastatest_treat_true + mastatest_treat_false),
mastatest_npv = mastatest_leave_true / (mastatest_leave_true + mastatest_leave_false),
) %>% select(biplate_treat_true:mastatest_npv)
ab_use_labels <- ab_use_intensive_stack %>% select(Biplate:mastatest_npv) %>%
pivot_longer(everything(),
names_to = "label",
values_to = "value")
ab_use_labels <- ab_use_labels %>% mutate(
index = case_when(
grepl('accumast',label,ignore.case = T) ~ "Accumast",
grepl('biplate',label,ignore.case = T) ~ "Biplate",
grepl('checkup',label,ignore.case = T) ~ "Check Up",
grepl('mastatest',label,ignore.case = T) ~ "Mastatest",
grepl('petrifilm',label,ignore.case = T) ~ "Petrifilm"
),
label_type = case_when(
grepl('ppv',label,ignore.case = T) ~ "ppv",
grepl('npv',label,ignore.case = T) ~ "npv",
TRUE ~ "total"
)
)
ab_use_long <- ab_use_intensive_stack %>%
select(biplate_treat_true:mastatest_leave_false) %>%
pivot_longer(everything(),names_to = "measure",values_to = "percent")
ab_use_long_intensive <- ab_use_long %>% mutate(
index = case_when(
grepl('accumast',measure) ~ "Accumast",
grepl('biplate',measure) ~ "Biplate",
grepl('checkup',measure) ~ "Check Up",
grepl('mastatest',measure) ~ "Mastatest",
grepl('petrifilm',measure) ~ "Petrifilm"
),
abx = case_when(
grepl('treat',measure) ~ "Antibiotic candidate",
grepl('leave',measure) ~ "Non-antibiotic candidate",
),
true = case_when(
grepl('true',measure) ~ "Correct classification",
grepl('false',measure) ~ "Incorrect classification",
)
)
ab_use_long_intensive_labels <- ab_use_long_intensive %>% group_by(index,abx) %>% summarise(
sum_allocated = sum(percent),
correct_allocated = ifelse(true=="Correct classification",percent,0),
percent_correct = round((correct_allocated/sum_allocated)*100,1)
) %>% filter(correct_allocated > 0)
# ab_use_long <- merge(ab_use_long,ab_use_labels,
# by="index")
ab_use_long_intensive <- merge(ab_use_long_intensive,
ab_use_long_intensive_labels,
by=c("index","abx"))This figure shows the percentage of cases identified as antibiotic treatment candidates (left panel) and non-antibiotic candidates (right) using 5 point of care testing systems. The green bars indicate the cases correctly classified and the red bars indicate cases incorrectly classified. For example, the Accumast testing system found that 23% of cases were antibiotic treatment candidates (i.e., Gram-positive pathogens). Of these, 66 (white number) were correctly identified as antibiotic treatment candidates.
ggplot() +
facet_wrap(~ abx, ncol = 7) +
geom_bar(aes(y = percent, x = index, fill = forcats::fct_rev(true)), data = ab_use_long_intensive, stat = "identity") +
geom_text(aes(y = sum_allocated, x = index, label = paste0(round(sum_allocated,0), "%")),
data = ab_use_long_intensive, vjust = -0.5, size = 4) +
geom_text(aes(y = 1, x = index, label = paste0(round(percent_correct,0), "%")),
color = "white", data = ab_use_long_intensive, vjust = -0.5, size = 4) +
theme(legend.position = "bottom",
legend.direction = "horizontal",
panel.background = element_blank(),
legend.title = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman", size = 12), # Adjust the overall text size
axis.ticks.x = element_blank(),
plot.title = element_text(size = 16), # Adjust the size of the plot title
axis.text = element_text(size = 10), # Adjust the size of the axis labels
axis.title = element_text(size = 12), # Adjust the size of the axis titles
strip.text = element_text(size = 14), # Adjust the size of the facet labels
legend.text = element_text(size = 10)) + # Adjust the size of the legend labels
ylim(0, 100) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20), minor_breaks = seq(0, 100, 5)) +
labs(x = "", y = "Percentage") +
ggtitle("Intensive herds")
ggsave("intensive_abx_use.tiff", width = 8, height = 5, dpi = 200)ab_use_non_intensive_stack <- data %>% filter(contamination==0 & system=="Non-intensive") %>% summarise(
culture_result = "Non-contaminated samples",
sample_n= sum(contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate_tx1_dz1 = sum(biplate_tx_gram_pos==1 & gram_pos==1,na.rm=T),
biplate_tx0_dz1 = sum(biplate_tx_gram_pos==0 & gram_pos==1,na.rm=T),
biplate_tx1_dz0 = sum(biplate_tx_gram_pos==1 & gram_pos==0,na.rm=T),
biplate_tx0_dz0 = sum(biplate_tx_gram_pos==0 & gram_pos==0,na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast_tx1_dz1 = sum(accumast_tx_gram_pos==1 & gram_pos==1,na.rm=T),
accumast_tx0_dz1 = sum(accumast_tx_gram_pos==0 & gram_pos==1,na.rm=T),
accumast_tx1_dz0 = sum(accumast_tx_gram_pos==1 & gram_pos==0,na.rm=T),
accumast_tx0_dz0 = sum(accumast_tx_gram_pos==0 & gram_pos==0,na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup_tx1_dz1 = sum(checkup_tx_gram_pos==1 & gram_pos==1,na.rm=T),
checkup_tx0_dz1 = sum(checkup_tx_gram_pos==0 & gram_pos==1,na.rm=T),
checkup_tx1_dz0 = sum(checkup_tx_gram_pos==1 & gram_pos==0,na.rm=T),
checkup_tx0_dz0 = sum(checkup_tx_gram_pos==0 & gram_pos==0,na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm_tx1_dz1 = sum(petrifilm_tx_gram_pos==1 & gram_pos==1,na.rm=T),
petrifilm_tx0_dz1 = sum(petrifilm_tx_gram_pos==0 & gram_pos==1,na.rm=T),
petrifilm_tx1_dz0 = sum(petrifilm_tx_gram_pos==1 & gram_pos==0,na.rm=T),
petrifilm_tx0_dz0 = sum(petrifilm_tx_gram_pos==0 & gram_pos==0,na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest_tx1_dz1 = sum(mastatest_tx_gram_pos==1 & gram_pos==1,na.rm=T),
mastatest_tx0_dz1 = sum(mastatest_tx_gram_pos==0 & gram_pos==1,na.rm=T),
mastatest_tx1_dz0 = sum(mastatest_tx_gram_pos==1 & gram_pos==0,na.rm=T),
mastatest_tx0_dz0 = sum(mastatest_tx_gram_pos==0 & gram_pos==0,na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
biplate_treat_true = round(100*biplate_tx1_dz1/(biplate_tx0+biplate_tx1),1),
biplate_treat_false = round(100*biplate_tx1_dz0/(biplate_tx0+biplate_tx1),1),
biplate_leave_true = round(100*biplate_tx0_dz0/(biplate_tx0+biplate_tx1),1),
biplate_leave_false = round(100*biplate_tx0_dz1/(biplate_tx0+biplate_tx1),1),
accumast_treat_true = round(100*accumast_tx1_dz1/(accumast_tx0+accumast_tx1),1),
accumast_treat_false = round(100*accumast_tx1_dz0/(accumast_tx0+accumast_tx1),1),
accumast_leave_true = round(100*accumast_tx0_dz0/(accumast_tx0+accumast_tx1),1),
accumast_leave_false = round(100*accumast_tx0_dz1/(accumast_tx0+accumast_tx1),1),
checkup_treat_true = round(100*checkup_tx1_dz1/(checkup_tx0+checkup_tx1),1),
checkup_treat_false = round(100*checkup_tx1_dz0/(checkup_tx0+checkup_tx1),1),
checkup_leave_true = round(100*checkup_tx0_dz0/(checkup_tx0+checkup_tx1),1),
checkup_leave_false = round(100*checkup_tx0_dz1/(checkup_tx0+checkup_tx1),1),
petrifilm_treat_true = round(100*petrifilm_tx1_dz1/(petrifilm_tx0+petrifilm_tx1),1),
petrifilm_treat_false = round(100*petrifilm_tx1_dz0/(petrifilm_tx0+petrifilm_tx1),1),
petrifilm_leave_true = round(100*petrifilm_tx0_dz0/(petrifilm_tx0+petrifilm_tx1),1),
petrifilm_leave_false = round(100*petrifilm_tx0_dz1/(petrifilm_tx0+petrifilm_tx1),1),
mastatest_treat_true = round(100*mastatest_tx1_dz1/(mastatest_tx0+mastatest_tx1),1),
mastatest_treat_false = round(100*mastatest_tx1_dz0/(mastatest_tx0+mastatest_tx1),1),
mastatest_leave_true = round(100*mastatest_tx0_dz0/(mastatest_tx0+mastatest_tx1),1),
mastatest_leave_false = round(100*mastatest_tx0_dz1/(mastatest_tx0+mastatest_tx1),1),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),
biplate_ppv = biplate_treat_true / (biplate_treat_true + biplate_treat_false),
biplate_npv = biplate_leave_true / (biplate_leave_true + biplate_leave_false),
accumast_ppv = accumast_treat_true / (accumast_treat_true + accumast_treat_false),
accumast_npv = accumast_leave_true / (accumast_leave_true + accumast_leave_false),
checkup_ppv = checkup_treat_true / (checkup_treat_true + checkup_treat_false),
checkup_npv = checkup_leave_true / (checkup_leave_true + checkup_leave_false),
petrifilm_ppv = petrifilm_treat_true / (petrifilm_treat_true + petrifilm_treat_false),
petrifilm_npv = petrifilm_leave_true / (petrifilm_leave_true + petrifilm_leave_false),
mastatest_ppv = mastatest_treat_true / (mastatest_treat_true + mastatest_treat_false),
mastatest_npv = mastatest_leave_true / (mastatest_leave_true + mastatest_leave_false),
) %>% select(biplate_treat_true:mastatest_npv)
ab_use_labels <- ab_use_non_intensive_stack %>% select(Biplate:mastatest_npv) %>%
pivot_longer(everything(),
names_to = "label",
values_to = "value")
ab_use_labels <- ab_use_labels %>% mutate(
index = case_when(
grepl('accumast',label,ignore.case = T) ~ "Accumast",
grepl('biplate',label,ignore.case = T) ~ "Biplate",
grepl('checkup',label,ignore.case = T) ~ "Check Up",
grepl('mastatest',label,ignore.case = T) ~ "Mastatest",
grepl('petrifilm',label,ignore.case = T) ~ "Petrifilm"
),
label_type = case_when(
grepl('ppv',label,ignore.case = T) ~ "ppv",
grepl('npv',label,ignore.case = T) ~ "npv",
TRUE ~ "total"
)
)
ab_use_long <- ab_use_non_intensive_stack %>%
select(biplate_treat_true:mastatest_leave_false) %>%
pivot_longer(everything(),names_to = "measure",values_to = "percent")
ab_use_long_non_intensive <- ab_use_long %>% mutate(
index = case_when(
grepl('accumast',measure) ~ "Accumast",
grepl('biplate',measure) ~ "Biplate",
grepl('checkup',measure) ~ "Check Up",
grepl('mastatest',measure) ~ "Mastatest",
grepl('petrifilm',measure) ~ "Petrifilm"
),
abx = case_when(
grepl('treat',measure) ~ "Antibiotic candidate",
grepl('leave',measure) ~ "Non-antibiotic candidate",
),
true = case_when(
grepl('true',measure) ~ "Correct classification",
grepl('false',measure) ~ "Incorrect classification",
)
)
ab_use_long_non_intensive_labels <- ab_use_long_non_intensive %>% group_by(index,abx) %>% summarise(
sum_allocated = sum(percent),
correct_allocated = ifelse(true=="Correct classification",percent,0),
percent_correct = round((correct_allocated/sum_allocated)*100,1)
) %>% filter(correct_allocated > 0)
# ab_use_long <- merge(ab_use_long,ab_use_labels,
# by="index")
ab_use_long_non_intensive <- merge(ab_use_long_non_intensive,
ab_use_long_non_intensive_labels,
by=c("index","abx"))This figure is the same as above, but for non-intensive herds.
ggplot() +
facet_wrap(~ abx, ncol = 7) +
geom_bar(aes(y = percent, x = index, fill = forcats::fct_rev(true)), data = ab_use_long_non_intensive, stat = "identity") +
geom_text(aes(y = sum_allocated, x = index, label = paste0(round(sum_allocated,0), "%")),
data = ab_use_long_non_intensive, vjust = -0.5, size = 4) +
geom_text(aes(y = 1, x = index, label = paste0(round(percent_correct,0), "%")),
color = "white", data = ab_use_long_non_intensive, vjust = -0.5, size = 4) +
theme(legend.position = "bottom",
legend.direction = "horizontal",
panel.background = element_blank(),
legend.title = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman", size = 12), # Adjust the overall text size
axis.ticks.x = element_blank(),
plot.title = element_text(size = 16), # Adjust the size of the plot title
axis.text = element_text(size = 10), # Adjust the size of the axis labels
axis.title = element_text(size = 12), # Adjust the size of the axis titles
strip.text = element_text(size = 14), # Adjust the size of the facet labels
legend.text = element_text(size = 10)) + # Adjust the size of the legend labels
ylim(0, 100) +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20), minor_breaks = seq(0, 100, 5)) +
labs(x = "", y = "Percentage") +
ggtitle("Non-intensive herds")
ggsave("non_intensive_abx_use.tiff", width = 8, height = 5, dpi = 200)tb_all_samples <- data %>% filter(contamination==0 | contamination==1) %>% summarise(
culture_result = "All samples",
sample_n= sum(contamination==1 | contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate_tested = biplate_tx1 + biplate_tx0,
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast_tested = accumast_tx1 + accumast_tx0,
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup_tested = checkup_tx1 + checkup_tx0,
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm_tested = petrifilm_tx1 + petrifilm_tx0,
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest_tested = mastatest_tx1 + mastatest_tx0,
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest,
biplate_tested,accumast_tested,checkup_tested,petrifilm_tested,mastatest_tested)
tb_all_samples_tested <- tb_all_samples %>% select(biplate_tested:mastatest_tested)
tb_all_samples <- tb_all_samples %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
tb_contaminated <- data %>% filter(contamination==1) %>% summarise(
culture_result = "Contaminated samples",
sample_n= sum(contamination==1, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
tb_non_contaminated <- data %>% filter(contamination==0) %>% summarise(
culture_result = "Non-contaminated samples",
sample_n= sum(contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
get_table <- function(path,pathname){
path <- enquo(path)
data %>% filter(contamination==0,!!path==1) %>% summarise(
culture_result = pathname,
sample_n= sum(!!path==1,na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate_tx_percent = (biplate_tx1 / (biplate_tx0 + biplate_tx1))*100,
accumast_tx1 = sum(accumast_tx_gram_pos==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos==0, na.rm=T),
accumast_tx_percent = (accumast_tx1 / (accumast_tx0 + accumast_tx1))*100,
checkup_tx1 = sum(checkup_tx_gram_pos==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos==0, na.rm=T),
checkup_tx_percent = (checkup_tx1 / (checkup_tx0 + checkup_tx1))*100,
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm_tx_percent = (petrifilm_tx1 / (petrifilm_tx0 + petrifilm_tx1))*100,
mastatest_tx1 = sum(mastatest_tx_gram_pos==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos==0, na.rm=T),
mastatest_tx_percent = (mastatest_tx1 / (mastatest_tx0 + mastatest_tx1))*100,
biplate = paste0(biplate_tx1, " (", round(biplate_tx_percent,1), "%)"),
accumast = paste0(accumast_tx1, " (", round(accumast_tx_percent,1), "%)"),
checkup = paste0(checkup_tx1, " (", round(checkup_tx_percent,1), "%)"),
petrifilm = paste0(petrifilm_tx1, " (", round(petrifilm_tx_percent,1), "%)"),
mastatest = paste0(mastatest_tx1, " (", round(mastatest_tx_percent,1), "%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result, sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
}
tb_gram_pos <- get_table(gram_pos,"Gram-positive")
tb_non_gram_pos <- get_table(non_gram_pos,"Non Gram-positive")
tb_gram_pos_maj <- get_table(gram_pos_maj,"Gram-positive major")
tb_sslo <- get_table(sslo,"Strep and Strep-like organisms")
tb_gram_pos_kleb <- get_table(gram_pos_kleb,"Gram-positive and Klebsiella")
tb_staph_aureus <- get_table(staph_aureus,"Staphylococcus aureus")
tb_strep_uberis <- get_table(strep_uberis,"Streptococcus uberis")
tb_strep_dys <- get_table(strep_dys,"Streptococcus dysgalactiae")
tb_gram_neg <- get_table(gram_neg,"Gram-negative")
tb_e_coli <- get_table(e_coli,"Escherichia coli")
tb_klebsiella <- get_table(klebsiella,"Klebsiella sp.")
tb_no_growth <- get_table(no_growth,"No growth")
table1 <- rbind(tb_all_samples,tb_contaminated,tb_non_contaminated,tb_gram_pos,tb_non_gram_pos,tb_gram_pos_maj,tb_gram_pos_kleb,tb_sslo,
tb_staph_aureus,tb_strep_uberis,tb_strep_dys,tb_gram_neg,tb_e_coli,tb_klebsiella,tb_no_growth)
klebsiella_not_targeted <- tb_klebsiella
write.csv(table1,"abx use pathogens.csv")The shows the total number of samples tested for each system. Note that the accumast testing system has missing values due to the plates being unavailable at times during the study.
tb_all_samples_tested| biplate_tested | accumast_tested | checkup_tested | petrifilm_tested | mastatest_tested |
|---|---|---|---|---|
| 641 | 550 | 639 | 640 | 634 |
The table below indicates what proportion of each sample result type (e.g. contamination, Strep uberis) that would be treated using each of the index testing systems. For example, of the 84.2% of Gram positive samples were identified as treatment targets by the Biplate testing system, which is slightly more than the Accumast system, which identified 75.9% of those samples as treatment targets.
table1 %>% select(culture_result:mastatest) %>% datatable(filter = "top")table1_long <- table1 %>% pivot_longer(cols = Biplate:Mastatest, names_to = "Test", values_to = "Percent")
table1_long <- table1_long %>% mutate(
culture_result_recode = case_when(
culture_result == "Escherichia coli" ~ "E. coli",
culture_result == "Staphylococcus aureus" ~ "Staph. aureus",
culture_result == "Streptococcus uberis" ~ "Strep. uberis",
culture_result == "Streptococcus dysgalactiae" ~ "Strep. dysgalactiae",
culture_result == "Klebsiella sp." ~ "Klebsiella spp.",
culture_result == "No growth" ~ "No growth",
)
)
outcomes_of_interest <- c("No growth",
"Escherichia coli",
"Klebsiella sp.",
"Staphylococcus aureus",
"Streptococcus uberis",
"Streptococcus dysgalactiae")
outcomes_level <- c("No growth",
"E. coli",
"Klebsiella spp.",
"Staph. aureus",
"Strep. uberis",
"Strep. dysgalactiae")
table1_long$Percent <- table1_long$Percent / 100
plot1 <- ggplot(table1_long %>% filter(culture_result %in% outcomes_of_interest)) +
aes(y = Percent, x = factor(culture_result_recode, level = outcomes_level)) +
geom_bar(aes(fill = Test), position = "dodge", stat="identity") +
scale_y_continuous(breaks=seq(0, 1, 0.2),labels = scales::percent, limits=c(0, 1)) +
ggtitle(paste0("Probability of antibiotic treatment")) +
labs(y = "Percentage of cows", x = "") +
scale_fill_brewer(palette="Set1") +
# geom_text(aes(x = Test, y = Percent), size=3.5, nudge_y = .02 ,colour='black') +
labs(y = "", x = "") +
#geom_text(aes(x = culture_result, y = Percent, fill = Test, label=Percent), size=3.5, nudge_y = .02 ,colour='black')
theme(panel.border = element_rect(colour = "black", fill=NA,size=1),
title=element_text(size=12,face="bold",family="Times"),
axis.text=element_text(size=12,family="Times"),
axis.title=element_text(size=11,face="bold",family="Times"),
axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.2, linetype = 'solid', colour = "grey"),
panel.grid.minor = element_line(size = 0.2, linetype = 'solid',colour = "grey"),
legend.text = element_text(colour="black",size=8,face="bold",family="Times"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())This figure shows the same thing as the table above.
plot1
This is another representation of the same data.
ggplot() +
facet_wrap(~ culture_result_recode, ncol = 3) +
geom_bar(aes(y = Percent, x = Test), data = table1_long %>% filter(culture_result %in% outcomes_of_interest),
stat = "identity", position = "dodge") +
geom_text(aes(y = Percent, x = Test, label = paste0(round(Percent*100,0), "%")),
data = table1_long %>% filter(culture_result %in% outcomes_of_interest), vjust = -0.5, size = 2) +
theme(legend.position = "bottom",
legend.direction = "horizontal",
panel.background = element_blank(),
legend.title = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman", size = 12),
axis.ticks.x = element_blank(),
plot.title = element_text(size = 16),
axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
strip.text = element_text(size = 14),
legend.text = element_text(size = 10)) +
scale_y_continuous(labels = scales::percent_format()) +
labs(x = "", y = "Percentage (%)")
ggsave("percent treated by pathogen.tiff", width = 10, height = 5, dpi = 200)# Repeat with Klebsiella in treated group
tb_all_samples <- data %>% filter(contamination==0 | contamination==1) %>% summarise(
culture_result = "All samples",
sample_n = sum(contamination==1 | contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos_kleb==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos_kleb==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos_kleb==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos_kleb==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos_kleb==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos_kleb==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
tb_contaminated <- data %>% filter(contamination==1) %>% summarise(
culture_result = "Contaminated samples",
sample_n= sum(contamination==1, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos_kleb==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos_kleb==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos_kleb==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos_kleb==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos_kleb==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos_kleb==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
tb_non_contaminated <- data %>% filter(contamination==0) %>% summarise(
culture_result = "Non-contaminated samples",
sample_n= sum(contamination==0, na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate = paste0(biplate_tx1, " (",round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),"%)"),
accumast_tx1 = sum(accumast_tx_gram_pos_kleb==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos_kleb==0, na.rm=T),
accumast = paste0(accumast_tx1, " (",round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),"%)"),
checkup_tx1 = sum(checkup_tx_gram_pos_kleb==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos_kleb==0, na.rm=T),
checkup = paste0(checkup_tx1, " (",round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),"%)"),
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm = paste0(petrifilm_tx1, " (",round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),"%)"),
mastatest_tx1 = sum(mastatest_tx_gram_pos_kleb==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos_kleb==0, na.rm=T),
mastatest = paste0(mastatest_tx1, " (",round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1),"%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result,sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
get_table <- function(path,pathname){
path <- enquo(path)
data %>% filter(contamination==0,!!path==1) %>% summarise(
culture_result = pathname,
sample_n= sum(!!path==1,na.rm=T),
biplate_tx1 = sum(biplate_tx_gram_pos==1, na.rm=T),
biplate_tx0 = sum(biplate_tx_gram_pos==0, na.rm=T),
biplate_tx_percent = (biplate_tx1 / (biplate_tx0 + biplate_tx1))*100,
accumast_tx1 = sum(accumast_tx_gram_pos_kleb==1, na.rm=T),
accumast_tx0 = sum(accumast_tx_gram_pos_kleb==0, na.rm=T),
accumast_tx_percent = (accumast_tx1 / (accumast_tx0 + accumast_tx1))*100,
checkup_tx1 = sum(checkup_tx_gram_pos_kleb==1, na.rm=T),
checkup_tx0 = sum(checkup_tx_gram_pos_kleb==0, na.rm=T),
checkup_tx_percent = (checkup_tx1 / (checkup_tx0 + checkup_tx1))*100,
petrifilm_tx1 = sum(petrifilm_tx_gram_pos==1, na.rm=T),
petrifilm_tx0 = sum(petrifilm_tx_gram_pos==0, na.rm=T),
petrifilm_tx_percent = (petrifilm_tx1 / (petrifilm_tx0 + petrifilm_tx1))*100,
mastatest_tx1 = sum(mastatest_tx_gram_pos_kleb==1, na.rm=T),
mastatest_tx0 = sum(mastatest_tx_gram_pos_kleb==0, na.rm=T),
mastatest_tx_percent = (mastatest_tx1 / (mastatest_tx0 + mastatest_tx1))*100,
biplate = paste0(biplate_tx1, " (", round(biplate_tx_percent,1), "%)"),
accumast = paste0(accumast_tx1, " (", round(accumast_tx_percent,1), "%)"),
checkup = paste0(checkup_tx1, " (", round(checkup_tx_percent,1), "%)"),
petrifilm = paste0(petrifilm_tx1, " (", round(petrifilm_tx_percent,1), "%)"),
mastatest = paste0(mastatest_tx1, " (", round(mastatest_tx_percent,1), "%)"),
Biplate = round(100*biplate_tx1/(biplate_tx0+biplate_tx1),1),
Accumast = round(100*accumast_tx1/(accumast_tx0+accumast_tx1),1),
Checkup = round(100*checkup_tx1/(checkup_tx0+checkup_tx1),1),
Petrifilm = round(100*petrifilm_tx1/(petrifilm_tx0+petrifilm_tx1),1),
Mastatest = round(100*mastatest_tx1/(mastatest_tx0+mastatest_tx1),1)
) %>% select(culture_result, sample_n,biplate,accumast,checkup,petrifilm,mastatest,
Biplate,Accumast,Checkup,Petrifilm,Mastatest)
}
tb_gram_pos <- get_table(gram_pos,"Gram-positive")
tb_non_gram_pos <- get_table(non_gram_pos,"Non Gram-positive")
tb_gram_pos_maj <- get_table(gram_pos_maj,"Gram-positive major")
tb_sslo <- get_table(sslo,"Strep and Strep-like organisms")
tb_gram_pos_kleb <- get_table(gram_pos_kleb,"Gram-positive and Klebsiella")
tb_staph_aureus <- get_table(staph_aureus,"Staphylococcus aureus")
tb_strep_uberis <- get_table(strep_uberis,"Streptococcus uberis")
tb_strep_dys <- get_table(strep_dys,"Streptococcus dysgalactiae")
tb_gram_neg <- get_table(gram_neg,"Gram-negative")
tb_e_coli <- get_table(e_coli,"Escherichia coli")
tb_klebsiella <- get_table(klebsiella,"Klebsiella sp.")
tb_no_growth <- get_table(no_growth,"No growth")
table2 <- rbind(tb_all_samples,tb_contaminated,tb_non_contaminated,tb_gram_pos,tb_non_gram_pos,tb_gram_pos_maj,tb_gram_pos_kleb,tb_sslo,
tb_staph_aureus,tb_strep_uberis,tb_strep_dys,tb_gram_neg,tb_e_coli,tb_klebsiella,tb_no_growth)
write.csv(table2,"abx use pathogens target includes klebsiella.csv")
klebsiella_targeted <- tb_klebsielladata %>% group_by(farm_number,system) %>% summarise(
n = n(),
klebsiella = sum(klebsiella==1,na.rm=T),
kleb_percent = klebsiella/n
) %>% filter(n > 9) %>% arrange(desc(kleb_percent))| farm_number | system | n | klebsiella | kleb_percent |
|---|---|---|---|---|
| 1 | Intensive | 47 | 20 | 0.4255319 |
| 2 | Intensive | 68 | 21 | 0.3088235 |
| 15 | Non-intensive | 19 | 2 | 0.1052632 |
| 19 | Intensive | 49 | 5 | 0.1020408 |
| 4 | Non-intensive | 13 | 1 | 0.0769231 |
| 24 | Intensive | 14 | 1 | 0.0714286 |
| 25 | Non-intensive | 26 | 1 | 0.0384615 |
| 20 | Non-intensive | 49 | 1 | 0.0204082 |
| 3 | Non-intensive | 15 | 0 | 0.0000000 |
| 5 | Non-intensive | 23 | 0 | 0.0000000 |
| 8 | Intensive | 11 | 0 | 0.0000000 |
| 9 | Non-intensive | 15 | 0 | 0.0000000 |
| 10 | Non-intensive | 82 | 0 | 0.0000000 |
| 11 | Non-intensive | 16 | 0 | 0.0000000 |
| 12 | Intensive | 18 | 0 | 0.0000000 |
| 13 | Non-intensive | 27 | 0 | 0.0000000 |
| 18 | Non-intensive | 19 | 0 | 0.0000000 |
| 22 | Non-intensive | 15 | 0 | 0.0000000 |
| 23 | Non-intensive | 16 | 0 | 0.0000000 |
| 26 | Non-intensive | 21 | 0 | 0.0000000 |
| 28 | Non-intensive | 33 | 0 | 0.0000000 |
klebsiella_targeted$method <- "Targeting Klebsiella"
klebsiella_not_targeted$method <- "Not Targeting Klebsiella"
kleb <- rbind(klebsiella_not_targeted,
klebsiella_targeted)
kleb %>% select(method,sample_n,accumast,checkup,mastatest)| method | sample_n | accumast | checkup | mastatest |
|---|---|---|---|---|
| Not Targeting Klebsiella | 52 | 7 (16.7%) | 19 (38%) | 9 (17.6%) |
| Targeting Klebsiella | 52 | 24 (57.1%) | 37 (71.2%) | 9 (17.6%) |
The table below indicates what proportion of each sample result type (e.g. contamination, Strep uberis) that would be treated using each of the index testing systems. For example, of the 71.6% of Gram positive or Klebsiella samples were identified as treatment targets by the Biplate testing system, which is slightly less than the Accumast system, which identified 73.7% of those samples as treatment targets.
table2 %>% select(culture_result:mastatest) %>% datatable(filter = "top")table2_long <- table2 %>% pivot_longer(cols = Biplate:Mastatest, names_to = "Test", values_to = "Percent")
table2_long <- table2_long %>% mutate(
culture_result_recode = case_when(
culture_result == "Escherichia coli" ~ "E. coli",
culture_result == "Staphylococcus aureus" ~ "Staph. aureus",
culture_result == "Streptococcus uberis" ~ "Strep. uberis",
culture_result == "Streptococcus dysgalactiae" ~ "Strep. dysgalactiae",
culture_result == "Klebsiella sp." ~ "Klebsiella spp.",
culture_result == "No growth" ~ "No growth",
)
)
outcomes_of_interest <- c("No growth",
"Escherichia coli",
"Klebsiella sp.",
"Staphylococcus aureus",
"Streptococcus uberis",
"Streptococcus dysgalactiae")
outcomes_level <- c("No growth",
"E. coli",
"Klebsiella spp.",
"Staph. aureus",
"Strep. uberis",
"Strep. dysgalactiae")
table2_long$Percent <- table2_long$Percent / 100This figure is the same data as the table above.
ggplot() +
facet_wrap(~ culture_result_recode, ncol = 3) +
geom_bar(aes(y = Percent, x = Test), data = table2_long %>% filter(culture_result %in% outcomes_of_interest),
stat = "identity", position = "dodge") +
geom_text(aes(y = Percent, x = Test, label = paste0(round(Percent*100,0), "%")),
data = table2_long %>% filter(culture_result %in% outcomes_of_interest), vjust = -0.5, size = 2) +
theme(legend.position = "bottom",
legend.direction = "horizontal",
panel.background = element_blank(),
legend.title = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman", size = 12),
axis.ticks.x = element_blank(),
plot.title = element_text(size = 16),
axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
strip.text = element_text(size = 14),
legend.text = element_text(size = 10)) +
scale_y_continuous(labels = scales::percent_format()) +
labs(x = "", y = "Percentage (%)")
ggsave("percent treated by pathogen when targeting klebsiella.tiff", width = 10, height = 5, dpi = 200)#Create a coliform group
rbind(data$iso1_dx %>% unique,
data$iso2_dx %>% unique) %>% as.vector %>% unique %>% write.csv('allpath.csv')
coliform_list <- c("Citrobacter freundii",
"Escherichia coli",
"Klebsiella oxytoca",
"Klebsiella pneumoniae",
"Klebsiella sp.",
"Klebsiella variicola",
"Serratia liquefaciens",
"Serratia marcescens")
data <- data %>% mutate(
coliform = case_when(
contamination== 0 & (iso1_dx %in% coliform_list | iso2_dx %in% coliform_list) ~ 1,
contamination== 0 ~ 0
),
accumast_coliform_aggregate = case_when(
accumast_e_coli == 1 |
accumast_enterobacter == 1 |
accumast_klebsiella == 1 ~ 1,
!is.na(accumast_contamination) ~ 0),
checkup_coliform_aggregate = case_when(
checkup_coliform_other == 1 |
checkup_e_coli == 1 |
checkup_klebsiella == 1 ~ 1,
!is.na(checkup_contamination) ~ 0
),
mastatest_coliform_aggregate = mastatest_coliform
)
# data %>% tabyl(coliform,contamination)
# data %>% tabyl(accumast_e_coli,accumast_coliform_aggregate)
# data %>% tabyl(accumast_klebsiella,accumast_coliform_aggregate)
# data %>% tabyl(accumast_enterobacter,accumast_coliform_aggregate)
# data %>% tabyl(coliform,accumast_coliform_aggregate)
# data %>% tabyl(coliform,checkup_coliform_aggregate)
# data %>% tabyl(coliform,mastatest_coliform_aggregate)test_eval <- function(reference,index){
reference <- enquo(reference)
index <- enquo(index)
dat <- data %>%
mutate(dis = factor(!!reference, levels = c(1,0), labels = c("Dis+","Dis-"))) %>%
mutate(tes = factor(!!index, levels = c(1,0), labels = c("Test+","Test-"))) %>%
group_by(tes, dis) %>%
summarise(n = n()) %>% drop_na
out <- epi.tests(dat, method = "exact", digits = 2,
conf.level = 0.95)
se <- out[["detail"]] %>% filter(statistic=="se") %>% select(est) %>% as.numeric %>% round(2)
sp <- out[["detail"]] %>% filter(statistic=="sp") %>% select(est) %>% as.numeric %>% round(2)
ppv <- out[["detail"]] %>% filter(statistic=="pv.pos") %>% select(est) %>% as.numeric %>% round(2)
npv <- out[["detail"]] %>% filter(statistic=="pv.neg") %>% select(est) %>% as.numeric %>% round(2)
target <- deparse(substitute(reference))
index <- deparse(substitute(index))
cbind(target,index,se,sp,ppv,npv)
}
test_eval_intensive <- function(reference,index){
reference <- enquo(reference)
index <- enquo(index)
dat <- data %>% filter(system=="Intensive") %>%
mutate(dis = factor(!!reference, levels = c(1,0), labels = c("Dis+","Dis-"))) %>%
mutate(tes = factor(!!index, levels = c(1,0), labels = c("Test+","Test-"))) %>%
group_by(tes, dis) %>%
summarise(n = n()) %>% drop_na
out <- epi.tests(dat, method = "exact", digits = 2,
conf.level = 0.95)
se <- out[["detail"]] %>% filter(statistic=="se") %>% select(est) %>% as.numeric %>% round(2)
sp <- out[["detail"]] %>% filter(statistic=="sp") %>% select(est) %>% as.numeric %>% round(2)
ppv <- out[["detail"]] %>% filter(statistic=="pv.pos") %>% select(est) %>% as.numeric %>% round(2)
npv <- out[["detail"]] %>% filter(statistic=="pv.neg") %>% select(est) %>% as.numeric %>% round(2)
target <- deparse(substitute(reference))
index <- deparse(substitute(index))
cbind(target,index,se,sp,ppv,npv)
}
test_eval_non_intensive <- function(reference,index){
reference <- enquo(reference)
index <- enquo(index)
dat <- data %>% filter(system=="Non-intensive") %>%
mutate(dis = factor(!!reference, levels = c(1,0), labels = c("Dis+","Dis-"))) %>%
mutate(tes = factor(!!index, levels = c(1,0), labels = c("Test+","Test-"))) %>%
group_by(tes, dis) %>%
summarise(n = n()) %>% drop_na
out <- epi.tests(dat, method = "exact", digits = 2,
conf.level = 0.95)
se <- out[["detail"]] %>% filter(statistic=="se") %>% select(est) %>% as.numeric %>% round(2)
sp <- out[["detail"]] %>% filter(statistic=="sp") %>% select(est) %>% as.numeric %>% round(2)
ppv <- out[["detail"]] %>% filter(statistic=="pv.pos") %>% select(est) %>% as.numeric %>% round(2)
npv <- out[["detail"]] %>% filter(statistic=="pv.neg") %>% select(est) %>% as.numeric %>% round(2)
target <- deparse(substitute(reference))
index <- deparse(substitute(index))
cbind(target,index,se,sp,ppv,npv)
}gram_pos_biplate <- test_eval(reference=gram_pos,index=biplate_gram_pos)
gram_pos_petrifilm <- test_eval(reference=gram_pos,index=petrifilm_gram_pos)
gram_pos_checkup <- test_eval(reference=gram_pos,index=checkup_gram_pos)
gram_pos_accumast <- test_eval(reference=gram_pos,index=accumast_gram_pos)
gram_pos_mastatest <- test_eval(reference=gram_pos,index=mastatest_gram_pos)
staph_aureus_checkup <- test_eval(reference=staph_aureus,index=checkup_staph_aureus)
staph_aureus_accumast <- test_eval(reference=staph_aureus,index=accumast_staph_aureus)
staph_aureus_mastatest <- test_eval(reference=staph_aureus,index=mastatest_staph_aureus)
strep_uberis_checkup <- test_eval(reference=strep_uberis,index=checkup_strep_uberis)
strep_uberis_accumast <- test_eval(reference=strep_uberis,index=accumast_strep_uberis)
strep_uberis_mastatest <- test_eval(reference=strep_uberis,index=mastatest_strep_uberis)
e_coli_checkup <- test_eval(reference=e_coli,index=checkup_e_coli)
e_coli_accumast <- test_eval(reference=e_coli,index=accumast_e_coli)
e_coli_mastatest <- test_eval(reference=e_coli,index=mastatest_e_coli)
klebsiella_checkup <- test_eval(reference=klebsiella,index=checkup_klebsiella)
klebsiella_accumast <- test_eval(reference=klebsiella,index=accumast_klebsiella)
klebsiella_mastatest <- test_eval(reference=klebsiella,index=mastatest_klebsiella)
coliform_checkup <- test_eval(reference=coliform,index=checkup_coliform_aggregate)
coliform_accumast <- test_eval(reference=coliform,index=accumast_coliform_aggregate)
coliform_mastatest <- test_eval(reference=coliform,index=mastatest_coliform_aggregate)
strep_dys_checkup <- test_eval(reference=strep_dys,index=checkup_strep_dys)
#strep_dys_accumast <- test_eval(reference=strep_dys,index=accumast_strep_dys)
strep_dys_mastatest <- test_eval(reference=strep_dys,index=mastatest_strep_dys)
table3 <- rbind(gram_pos_checkup,gram_pos_accumast,gram_pos_mastatest,
coliform_checkup,coliform_accumast,coliform_mastatest,
staph_aureus_checkup,staph_aureus_accumast,staph_aureus_mastatest,
strep_uberis_checkup,strep_uberis_accumast,strep_uberis_mastatest,
e_coli_checkup,e_coli_accumast,
klebsiella_checkup,klebsiella_accumast,
strep_dys_checkup) %>% as.data.frame()
table3 <- table3 %>% mutate(
index_test = case_when(
str_detect(index, "checkup") ~ "Checkup",
str_detect(index, "accumast") ~ "Accumast",
str_detect(index, "mastatest") ~ "Mastatest",
)
)
table3 <- table3 %>% select(target,index_test,se,sp,ppv,npv)# table3 %>%
# pivot_wider(names_from = index_test, values_from = c(se,sp)) %>%
# datatable(filter = "top")gram_pos_biplate <- test_eval_non_intensive(reference=gram_pos,index=biplate_gram_pos)
gram_pos_petrifilm <- test_eval_non_intensive(reference=gram_pos,index=petrifilm_gram_pos)
gram_pos_checkup <- test_eval_non_intensive(reference=gram_pos,index=checkup_gram_pos)
gram_pos_accumast <- test_eval_non_intensive(reference=gram_pos,index=accumast_gram_pos)
gram_pos_mastatest <- test_eval_non_intensive(reference=gram_pos,index=mastatest_gram_pos)
staph_aureus_checkup <- test_eval_non_intensive(reference=staph_aureus,index=checkup_staph_aureus)
staph_aureus_accumast <- test_eval_non_intensive(reference=staph_aureus,index=accumast_staph_aureus)
staph_aureus_mastatest <- test_eval_non_intensive(reference=staph_aureus,index=mastatest_staph_aureus)
strep_uberis_checkup <- test_eval_non_intensive(reference=strep_uberis,index=checkup_strep_uberis)
strep_uberis_accumast <- test_eval_non_intensive(reference=strep_uberis,index=accumast_strep_uberis)
strep_uberis_mastatest <- test_eval_non_intensive(reference=strep_uberis,index=mastatest_strep_uberis)
e_coli_checkup <- test_eval_non_intensive(reference=e_coli,index=checkup_e_coli)
e_coli_accumast <- test_eval_non_intensive(reference=e_coli,index=accumast_e_coli)
#e_coli_mastatest <- test_eval_non_intensive(reference=e_coli,index=mastatest_e_coli)
klebsiella_checkup <- test_eval_non_intensive(reference=klebsiella,index=checkup_klebsiella)
klebsiella_accumast <- test_eval_non_intensive(reference=klebsiella,index=accumast_klebsiella)
#klebsiella_mastatest <- test_eval_non_intensive(reference=klebsiella,index=mastatest_klebsiella)
coliform_checkup <- test_eval_non_intensive(reference=coliform,index=checkup_coliform_aggregate)
coliform_accumast <- test_eval_non_intensive(reference=coliform,index=accumast_coliform_aggregate)
coliform_mastatest <- test_eval_non_intensive(reference=coliform,index=mastatest_coliform_aggregate)
strep_dys_checkup <- test_eval_non_intensive(reference=strep_dys,index=checkup_strep_dys)
#strep_dys_accumast <- test_eval_non_intensive(reference=strep_dys,index=accumast_strep_dys)
#strep_dys_mastatest <- test_eval_non_intensive(reference=strep_dys,index=mastatest_strep_dys)
table4 <- rbind(gram_pos_checkup,gram_pos_accumast,gram_pos_mastatest,
coliform_checkup,coliform_accumast,coliform_mastatest,
staph_aureus_checkup,staph_aureus_accumast,staph_aureus_mastatest,
strep_uberis_checkup,strep_uberis_accumast,strep_uberis_mastatest,
e_coli_checkup,e_coli_accumast,
klebsiella_checkup,klebsiella_accumast,
strep_dys_checkup) %>% as.data.frame()
table4 <- table4 %>% mutate(
index_test = case_when(
str_detect(index, "checkup") ~ "Checkup",
str_detect(index, "accumast") ~ "Accumast",
str_detect(index, "mastatest") ~ "Mastatest",
)
)
table4$ppv_non_intensive <- table4$ppv
table4$npv_non_intensive <- table4$npv
table4 <- table4 %>% select(target,index_test,ppv_non_intensive,npv_non_intensive)gram_pos_biplate <- test_eval_intensive(reference=gram_pos,index=biplate_gram_pos)
gram_pos_petrifilm <- test_eval_intensive(reference=gram_pos,index=petrifilm_gram_pos)
gram_pos_checkup <- test_eval_intensive(reference=gram_pos,index=checkup_gram_pos)
gram_pos_accumast <- test_eval_intensive(reference=gram_pos,index=accumast_gram_pos)
gram_pos_mastatest <- test_eval_intensive(reference=gram_pos,index=mastatest_gram_pos)
staph_aureus_checkup <- test_eval_intensive(reference=staph_aureus,index=checkup_staph_aureus)
staph_aureus_accumast <- test_eval_intensive(reference=staph_aureus,index=accumast_staph_aureus)
staph_aureus_mastatest <- test_eval_intensive(reference=staph_aureus,index=mastatest_staph_aureus)
strep_uberis_checkup <- test_eval_intensive(reference=strep_uberis,index=checkup_strep_uberis)
strep_uberis_accumast <- test_eval_intensive(reference=strep_uberis,index=accumast_strep_uberis)
strep_uberis_mastatest <- test_eval_intensive(reference=strep_uberis,index=mastatest_strep_uberis)
e_coli_checkup <- test_eval_intensive(reference=e_coli,index=checkup_e_coli)
e_coli_accumast <- test_eval_intensive(reference=e_coli,index=accumast_e_coli)
#e_coli_mastatest <- test_eval_intensive(reference=e_coli,index=mastatest_e_coli)
klebsiella_checkup <- test_eval_intensive(reference=klebsiella,index=checkup_klebsiella)
klebsiella_accumast <- test_eval_intensive(reference=klebsiella,index=accumast_klebsiella)
#klebsiella_mastatest <- test_eval_intensive(reference=klebsiella,index=mastatest_klebsiella)
coliform_checkup <- test_eval_intensive(reference=coliform,index=checkup_coliform_aggregate)
coliform_accumast <- test_eval_intensive(reference=coliform,index=accumast_coliform_aggregate)
coliform_mastatest <- test_eval_intensive(reference=coliform,index=mastatest_coliform_aggregate)
strep_dys_checkup <- test_eval_intensive(reference=strep_dys,index=checkup_strep_dys)
#strep_dys_accumast <- test_eval_intensive(reference=strep_dys,index=accumast_strep_dys)
#strep_dys_mastatest <- test_eval_intensive(reference=strep_dys,index=mastatest_strep_dys)
table5 <- rbind(gram_pos_checkup,gram_pos_accumast,gram_pos_mastatest,
coliform_checkup,coliform_accumast,coliform_mastatest,
staph_aureus_checkup,staph_aureus_accumast,staph_aureus_mastatest,
strep_uberis_checkup,strep_uberis_accumast,strep_uberis_mastatest,
e_coli_checkup,e_coli_accumast,
klebsiella_checkup,klebsiella_accumast,
strep_dys_checkup) %>% as.data.frame()
table5 <- table5 %>% mutate(
index_test = case_when(
str_detect(index, "checkup") ~ "Checkup",
str_detect(index, "accumast") ~ "Accumast",
str_detect(index, "mastatest") ~ "Mastatest",
)
)
table5$ppv_intensive <- table5$ppv
table5$npv_intensive <- table5$npv
table5 <- table5 %>% select(target,index_test,ppv_intensive,npv_intensive)
test_characteristics <- merge(table3,table4,
by=c("target","index_test"))
test_characteristics <- merge(test_characteristics,table5,
by=c("target","index_test"))The table below shows the sensitivity and specificity for selected index tests when attempting to identify Gram-positive pathogens, E. coli, Klebsiella spp., Staph aureus, Strep dysgalactiae, and Strep uberis. The table also outlines the predictive values in intensive and non-intensive herds, which have significantly different pathogen profiles.
prev_table <- data %>% group_by(system) %>% summarise(
coliform = sum(coliform==1,na.rm=T)/sum(!is.na(coliform),na.rm=T),
e_coli = sum(e_coli==1,na.rm=T)/sum(!is.na(e_coli),na.rm=T),
gram_pos = sum(gram_pos==1,na.rm=T)/sum(!is.na(gram_pos),na.rm=T),
klebsiella = sum(klebsiella==1,na.rm=T)/sum(!is.na(klebsiella),na.rm=T),
staph_aureus = sum(staph_aureus==1,na.rm=T)/sum(!is.na(staph_aureus),na.rm=T),
strep_dys = sum(strep_dys==1,na.rm=T)/sum(!is.na(strep_dys),na.rm=T),
strep_uberis = sum(strep_uberis==1,na.rm=T)/sum(!is.na(strep_uberis),na.rm=T),
) %>% filter(!is.na(system))
prev_table <- prev_table %>% pivot_longer(coliform:strep_uberis,
names_to = "target",
values_to = "prevalence"
)
prev_table <- prev_table %>% pivot_wider(names_from = system,
values_from = prevalence)
colnames(prev_table) <- c("target","prev_intensive","prev_non_intensive")
test_characteristics <- test_characteristics %>% mutate(
target = case_when(
target == "~coliform" ~ "coliform",
target == "~e_coli" ~ "e_coli",
target == "~gram_pos" ~ "gram_pos",
target == "~klebsiella" ~ "klebsiella",
target == "~staph_aureus" ~ "staph_aureus",
target == "~strep_dys" ~ "strep_dys",
target == "~strep_uberis" ~ "strep_uberis"
)
)
test_characteristics <- merge(test_characteristics %>% select(target,index_test,se,sp,ppv,npv),
prev_table,
by = "target")
test_characteristics <- test_characteristics %>% mutate(
se = as.numeric(se),
sp = as.numeric(sp),
ppv_intensive = round((se*prev_intensive) / ((se*prev_intensive) + (1-sp)*(1-prev_intensive)),2),
npv_intensive = round((sp*(1-prev_intensive)) / ((1-se)*prev_intensive + sp*(1-prev_intensive)),2),
ppv_non_intensive = round((se*prev_non_intensive) / ((se*prev_non_intensive) + (1-sp)*(1-prev_non_intensive)),2),
npv_non_intensive = round((sp*(1-prev_non_intensive)) / ((1-se)*prev_non_intensive + sp*(1-prev_non_intensive)),2)
) %>% select(-prev_intensive,-prev_non_intensive)write.csv(test_characteristics,"test characteristics.csv")
test_characteristics %>% datatable(filter = "top")The tables below show the what each index test found for samples identified by lab culture as: contaminated, no significant growth, E. coli, Klebsiella, Staph aureus, Strep dysgalactiae, and Strep uberis
data$all <- 1
summary <- data %>% select(contains("accumast"),strep_uberis)
vars <- colnames(summary) %>% as_tibble
#data %>% tabyl(accumast_contamination)
accumast_eval <- function(path,pathname){
path <- enquo(path)
data %>% filter(!!path==1) %>% summarise(
accumast_valid_n = sum(!is.na(accumast_contamination)==T),
cases = accumast_valid_n %>% as.character,
accumast_contamination = paste0(sum(accumast_contamination, na.rm=T),
" (",
round((sum(accumast_contamination, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_e_coli = paste0(sum(accumast_e_coli, na.rm=T),
" (",
round((sum(accumast_e_coli, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_enterobacter = paste0(sum(accumast_enterobacter, na.rm=T),
" (",
round((sum(accumast_enterobacter, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_enterococcus = paste0(sum(accumast_enterococcus, na.rm=T),
" (",
round((sum(accumast_enterococcus, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_gram_neg_non_typeable = paste0(sum(accumast_gram_neg_non_typeable, na.rm=T),
" (",
round((sum(accumast_gram_neg_non_typeable, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_klebsiella = paste0(sum(accumast_klebsiella, na.rm=T),
" (",
round((sum(accumast_klebsiella, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_lactococcus = paste0(sum(accumast_lactococcus, na.rm=T),
" (",
round((sum(accumast_lactococcus, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_no_growth = paste0(sum(accumast_no_growth, na.rm=T),
" (",
round((sum(accumast_no_growth, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_pseudomonas = paste0(sum(accumast_pseudomonas, na.rm=T),
" (",
round((sum(accumast_pseudomonas, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_sslo_non_typeable = paste0(sum(accumast_sslo_non_typeable, na.rm=T),
" (",
round((sum(accumast_sslo_non_typeable, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_staph_aureus = paste0(sum(accumast_staph_aureus, na.rm=T),
" (",
round((sum(accumast_staph_aureus, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_staph_chromogenes = paste0(sum(accumast_staph_chromogenes, na.rm=T),
" (",
round((sum(accumast_staph_chromogenes, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_staph_haemolyticus = paste0(sum(accumast_staph_haemolyticus, na.rm=T),
" (",
round((sum(accumast_staph_haemolyticus, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_staph_non_typeable = paste0(sum(accumast_staph_non_typeable, na.rm=T),
" (",
round((sum(accumast_staph_non_typeable, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_strep_spp = paste0(sum(accumast_strep_spp, na.rm=T),
" (",
round((sum(accumast_strep_spp, na.rm=T)/accumast_valid_n)*100,1),
"%)"),
accumast_strep_uberis = paste0(sum(accumast_strep_uberis, na.rm=T),
" (",
round((sum(accumast_strep_uberis, na.rm=T)/accumast_valid_n)*100,1),
"%)")
) %>% select(-accumast_valid_n) %>% pivot_longer(cols = everything(),
names_to = "Target",
values_to = pathname)
}
accumast_all <- accumast_eval(path=all,"All samples")
accumast_contamination <- accumast_eval(path=contamination,"Contamination")
accumast_no_growth <- accumast_eval(path=no_growth,"No growth")
accumast_e_coli <- accumast_eval(path=e_coli,"E. coli")
accumast_klebsiella <- accumast_eval(path=klebsiella,"Klebsiella spp.")
accumast_staph_aureus <- accumast_eval(path=staph_aureus,"Staph aureus")
accumast_strep_uberis <- accumast_eval(path=strep_uberis,"Strep uberis")
accumast_strep_dys <- accumast_eval(path=strep_dys,"Strep dys")
df_list <- list(
accumast_all,
accumast_contamination,
accumast_no_growth,
accumast_e_coli,
accumast_klebsiella,
accumast_staph_aureus,
accumast_strep_uberis,
accumast_strep_dys
)
accumast_out <- df_list %>% reduce(full_join, by='Target')
write.csv(accumast_out,'accumast out.csv')
accumast_out %>% datatable(filter = "top")#summarytools::dfSummary(summary %>% filter(strep_uberis==1),valid.col=FALSE, graph.magnif=0.8, style="grid")summary <- data %>% select(contains("biplate"),strep_uberis)
vars <- colnames(summary) %>% as_tibble
biplate_eval <- function(path,pathname){
path <- enquo(path)
data %>% filter(!!path==1) %>% summarise(
biplate_valid_n = sum(!is.na(biplate_contaminated)==T),
cases = biplate_valid_n %>% as.character,
biplate_contaminated = paste0(sum(biplate_contaminated, na.rm=T),
" (",
round((sum(biplate_contaminated, na.rm=T)/biplate_valid_n)*100,1),
"%)"),
biplate_no_growth = paste0(sum(biplate_no_growth, na.rm=T),
" (",
round((sum(biplate_no_growth, na.rm=T)/biplate_valid_n)*100,1),
"%)"),
biplate_gram_pos = paste0(sum(biplate_gram_pos, na.rm=T),
" (",
round((sum(biplate_gram_pos, na.rm=T)/biplate_valid_n)*100,1),
"%)"),
biplate_gram_neg = paste0(sum(biplate_gram_neg, na.rm=T),
" (",
round((sum(biplate_gram_neg, na.rm=T)/biplate_valid_n)*100,1),
"%)"),
) %>% select(-biplate_valid_n) %>% pivot_longer(cols = everything(),
names_to = "Target",
values_to = pathname)
}
biplate_all <- biplate_eval(path=all,"All samples")
biplate_contamination <- biplate_eval(path=contamination,"Contamination")
biplate_no_growth <- biplate_eval(path=no_growth,"No growth")
biplate_e_coli <- biplate_eval(path=e_coli,"E. coli")
biplate_klebsiella <- biplate_eval(path=klebsiella,"Klebsiella spp.")
biplate_staph_aureus <- biplate_eval(path=staph_aureus,"Staph aureus")
biplate_strep_uberis <- biplate_eval(path=strep_uberis,"Strep uberis")
biplate_strep_dys <- biplate_eval(path=strep_dys,"Strep dys")
df_list <- list(
biplate_all,
biplate_contamination,
biplate_no_growth,
biplate_e_coli,
biplate_klebsiella,
biplate_staph_aureus,
biplate_strep_uberis,
biplate_strep_dys
)
biplate_out <- df_list %>% reduce(full_join, by='Target')
write.csv(biplate_out,'biplate out.csv')
biplate_out %>% datatable(filter = "top")#summarytools::dfSummary(summary %>% filter(strep_uberis==1),valid.col=FALSE, graph.magnif=0.8, style="grid")data <- data %>% mutate(
checkup_no_growth = case_when(
checkup_aerobic_growth == 1 ~ 0,
checkup_aerobic_growth == 0 ~ 1
)
)
summary <- data %>% select(contains("checkup"),strep_uberis)
vars <- colnames(summary) %>% as_tibble
checkup_eval <- function(path,pathname){
path <- enquo(path)
data %>% filter(!!path==1) %>% summarise(
checkup_valid_n = sum(!is.na(checkup_contamination)==T),
cases = checkup_valid_n %>% as.character,
checkup_contaminated = paste0(sum(checkup_contamination, na.rm=T),
" (",
round((sum(checkup_contamination, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_coliform_other = paste0(sum(checkup_coliform_other, na.rm=T),
" (",
round((sum(checkup_coliform_other, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_e_coli = paste0(sum(checkup_e_coli, na.rm=T),
" (",
round((sum(checkup_e_coli, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_aerobic_non_typeable = paste0(sum(checkup_aerobic_non_typeable, na.rm=T),
" (",
round((sum(checkup_aerobic_non_typeable, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_gram_neg_non_typeable = paste0(sum(checkup_gram_neg_non_typeable, na.rm=T),
" (",
round((sum(checkup_gram_neg_non_typeable, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_gram_pos_non_typeable = paste0(sum(checkup_gram_pos_non_typeable, na.rm=T),
" (",
round((sum(checkup_gram_pos_non_typeable, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_klebsiella = paste0(sum(checkup_klebsiella, na.rm=T),
" (",
round((sum(checkup_klebsiella, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_no_growth = paste0(sum(checkup_no_growth, na.rm=T),
" (",
round((sum(checkup_no_growth, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_nas = paste0(sum(checkup_nas, na.rm=T),
" (",
round((sum(checkup_nas, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_pseudomonas = paste0(sum(checkup_pseudomonas, na.rm=T),
" (",
round((sum(checkup_pseudomonas, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_slo = paste0(sum(checkup_slo, na.rm=T),
" (",
round((sum(checkup_slo, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_staph_agressive = paste0(sum(checkup_staph_agressive, na.rm=T),
" (",
round((sum(checkup_staph_agressive, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_staph_aureus = paste0(sum(checkup_staph_aureus, na.rm=T),
" (",
round((sum(checkup_staph_aureus, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_strep_ag = paste0(sum(checkup_strep_ag, na.rm=T),
" (",
round((sum(checkup_strep_ag, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_strep_dys = paste0(sum(checkup_strep_dys, na.rm=T),
" (",
round((sum(checkup_strep_dys, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_strep_uberis = paste0(sum(checkup_strep_uberis, na.rm=T),
" (",
round((sum(checkup_strep_uberis, na.rm=T)/checkup_valid_n)*100,1),
"%)"),
checkup_yeast = paste0(sum(checkup_yeast, na.rm=T),
" (",
round((sum(checkup_yeast, na.rm=T)/checkup_valid_n)*100,1),
"%)")
) %>% select(-checkup_valid_n) %>% pivot_longer(cols = everything(),
names_to = "Target",
values_to = pathname)
}
checkup_all <- checkup_eval(path=all,"All samples")
checkup_contamination <- checkup_eval(path=contamination,"Contamination")
checkup_no_growth <- checkup_eval(path=no_growth,"No growth")
checkup_e_coli <- checkup_eval(path=e_coli,"E. coli")
checkup_klebsiella <- checkup_eval(path=klebsiella,"Klebsiella spp.")
checkup_staph_aureus <- checkup_eval(path=staph_aureus,"Staph aureus")
checkup_strep_uberis <- checkup_eval(path=strep_uberis,"Strep uberis")
checkup_strep_dys <- checkup_eval(path=strep_dys,"Strep dys")
df_list <- list(
checkup_all,
checkup_contamination,
checkup_no_growth,
checkup_e_coli,
checkup_klebsiella,
checkup_staph_aureus,
checkup_strep_uberis,
checkup_strep_dys
)
checkup_out <- df_list %>% reduce(full_join, by='Target')
write.csv(checkup_out,'checkup out.csv')
checkup_out %>% datatable(filter = "top")#summarytools::dfSummary(summary %>% filter(strep_uberis==1),valid.col=FALSE, graph.magnif=0.8, style="grid")
#check <- data %>% filter(no_growth==1) %>% select(record_id,checkup_contamination:checkup_no_growth)
#write.csv(check,'check.csv')mastatest_summary <- data %>% select(contains("mastatest"))
vars <- colnames(data) %>% as_tibble
mastatest_eval <- function(path,pathname){
path <- enquo(path)
data %>% filter(!!path==1) %>% summarise(
mastatest_valid_n = sum(mastatest_complete==2),
cases = mastatest_valid_n %>% as.character,
mastatest_nas = paste0(sum(mastatest_nas==1, na.rm=T),
" (",
round((sum(mastatest_nas==1, na.rm=T)/mastatest_valid_n)*100,1),
"%)"),
mastatest_coliform = paste0(sum(mastatest_coliform==1, na.rm=T),
" (",
round((sum(mastatest_coliform==1, na.rm=T)/mastatest_valid_n)*100,1),
"%)"),
mastatest_other_gram_pos = paste0(sum(mastatest_other_gram_pos==1, na.rm=T),
" (",
round((sum(mastatest_other_gram_pos==1, na.rm=T)/mastatest_valid_n)*100,1),
"%)"),
mastatest_strep_uberis = paste0(sum(mastatest_strep_uberis==1, na.rm=T),
" (",
round((sum(mastatest_strep_uberis==1, na.rm=T)/mastatest_valid_n)*100,1),
"%)"),
mastatest_staph_aureus = paste0(sum(mastatest_staph_aureus==1, na.rm=T),
" (",
round((sum(mastatest_staph_aureus==1, na.rm=T)/mastatest_valid_n)*100,1),
"%)"),
mastatest_other = paste0(sum(mastatest_other==1, na.rm=T),
" (",
round((sum(mastatest_other==1, na.rm=T)/mastatest_valid_n)*100,1),
"%)"),
mastatest_no_growth = paste0(sum(mastatest_no_growth, na.rm=T),
" (",
round((sum(mastatest_no_growth, na.rm=T)/mastatest_valid_n)*100,1),
"%)")
) %>% select(-mastatest_valid_n) %>% pivot_longer(cols = everything(),
names_to = "Target",
values_to = pathname)
}
mastatest_all <- mastatest_eval(path=all,"All samples")
mastatest_contamination <- mastatest_eval(path=contamination,"Contamination")
mastatest_no_growth <- mastatest_eval(path=no_growth,"No growth")
mastatest_e_coli <- mastatest_eval(path=e_coli,"E. coli")
mastatest_klebsiella <- mastatest_eval(path=klebsiella,"Klebsiella spp.")
mastatest_staph_aureus <- mastatest_eval(path=staph_aureus,"Staph aureus")
mastatest_strep_uberis <- mastatest_eval(path=strep_uberis,"Strep uberis")
mastatest_strep_dys <- mastatest_eval(path=strep_dys,"Strep dys")
mastatest_strep_dys <- mastatest_eval(path=strep_dys,"Strep dys")
df_list <- list(
mastatest_all,
mastatest_contamination,
mastatest_no_growth,
mastatest_e_coli,
mastatest_klebsiella,
mastatest_staph_aureus,
mastatest_strep_uberis,
mastatest_strep_dys
)
mastatest_out <- df_list %>% reduce(full_join, by='Target')
write.csv(mastatest_out,'mastatest out.csv')
mastatest_out %>% datatable(filter = "top")#e_coli <- data %>% filter(e_coli==1) %>% select(contains("mastatest")) %>% filter()summary <- data %>% select(contains("petrifilm"),strep_uberis)
vars <- colnames(summary) %>% as_tibble
petrifilm_eval <- function(path,pathname){
path <- enquo(path)
data %>% filter(!!path==1) %>% summarise(
petrifilm_valid_n = sum(!is.na(petrifilm_no_growth)==T),
cases = petrifilm_valid_n %>% as.character,
petrifilm_no_growth = paste0(sum(petrifilm_no_growth, na.rm=T),
" (",
round((sum(petrifilm_no_growth, na.rm=T)/petrifilm_valid_n)*100,1),
"%)"),
petrifilm_gram_pos = paste0(sum(petrifilm_gram_pos, na.rm=T),
" (",
round((sum(petrifilm_gram_pos, na.rm=T)/petrifilm_valid_n)*100,1),
"%)"),
petrifilm_gram_neg = paste0(sum(petrifilm_gram_neg, na.rm=T),
" (",
round((sum(petrifilm_gram_neg, na.rm=T)/petrifilm_valid_n)*100,1),
"%)"),
petrifilm_invalid = paste0(sum(petrifilm_invalid_result, na.rm=T),
" (",
round((sum(petrifilm_invalid_result, na.rm=T)/petrifilm_valid_n)*100,1),
"%)"),
) %>% select(-petrifilm_valid_n) %>% pivot_longer(cols = everything(),
names_to = "Target",
values_to = pathname)
}
petrifilm_all <- petrifilm_eval(path=all,"All samples")
petrifilm_contamination <- petrifilm_eval(path=contamination,"Contamination")
petrifilm_no_growth <- petrifilm_eval(path=no_growth,"No growth")
petrifilm_e_coli <- petrifilm_eval(path=e_coli,"E. coli")
petrifilm_klebsiella <- petrifilm_eval(path=klebsiella,"Klebsiella spp.")
petrifilm_staph_aureus <- petrifilm_eval(path=staph_aureus,"Staph aureus")
petrifilm_strep_uberis <- petrifilm_eval(path=strep_uberis,"Strep uberis")
petrifilm_strep_dys <- petrifilm_eval(path=strep_dys,"Strep dys")
df_list <- list(
petrifilm_all,
petrifilm_contamination,
petrifilm_no_growth,
petrifilm_e_coli,
petrifilm_klebsiella,
petrifilm_staph_aureus,
petrifilm_strep_uberis,
petrifilm_strep_dys
)
petrifilm_out <- df_list %>% reduce(full_join, by='Target')
write.csv(petrifilm_out,'petrifilm out.csv')
petrifilm_out %>% datatable(filter = "top")#summarytools::dfSummary(summary %>% filter(strep_uberis==1),valid.col=FALSE, graph.magnif=0.8, style="grid")out <- rbind(accumast_out,
biplate_out,
checkup_out,
mastatest_out,
petrifilm_out)
out$`All samples` <- str_sub(out$`All samples`,4,-3) %>% parse_number / 100
out$Contamination <- str_sub(out$Contamination,4,-3) %>% parse_number / 100
out$`No growth` <- str_sub(out$`No growth`,4,-3) %>% parse_number / 100
out$`E. coli` <- str_sub(out$`E. coli`,4,-3) %>% parse_number / 100
out$`Klebsiella spp.` <- str_sub(out$`Klebsiella spp.`,4,-3) %>% parse_number / 100
out$`Staph aureus` <- str_sub(out$`Staph aureus`,4,-3) %>% parse_number / 100
out$`Strep uberis` <- str_sub(out$`Strep uberis`,4,-3) %>% parse_number / 100
out$`Strep dys` <- str_sub(out$`Strep dys`,4,-3) %>% parse_number / 100
out <- out %>% filter(Target != "cases")
out <- out %>% mutate(
index = case_when(
grepl('accumast',Target) ~ "Accumast",
grepl('biplate',Target) ~ "Biplate",
grepl('checkup',Target) ~ "Check Up",
grepl('mastatest',Target) ~ "Mastatest",
grepl('petrifilm',Target) ~ "Petrifilm"
)
)
out$target <- str_replace(out$Target, "^.*?_", "")
out$target <- str_replace_all(out$target, "_", " ")
out$target <- str_to_sentence(out$target)
out$index_result <- out$target
out <- out %>% select(index,index_result,Contamination:`Strep dys`)
#out %>% datatable(filter = "top")The plots below are another way of looking at the data shown above. Each bar chart shows the top 5 results for samples identified by lab culture as either: contaminated, no significant growth, E. coli, Klebseilla, Staph aureus, Strep dysgalactiae, and Strep uberis
tab <- out %>% pivot_longer(
cols = Contamination:`Strep dys`, values_to = "Percent", names_to = "Target"
)
# tab %>% bar_chart(y = value, x = Result, top_n = 3,facet = name)
library(cowplot)
top5 <- function(target,index){
target <- enquo(target)
index <- enquo(index)
tab %>%
filter(Target == !!target & index == !!index) %>%
bar_chart(y = Percent, x = index_result, top_n = 5) +
ggtitle(paste0(quo_name(index),
" results when lab culture finds: ",
quo_name(target)
)) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())
}
all_targets <- unique(tab$Target)
all_index <- unique(tab$index)
#top5("Contamination","Accumast")
#top3("Contamination","Accumast")
combinations <- expand.grid(target = all_targets, index = all_index) %>% arrange(target)
# apply the top4 function to each combination of target and index# tab2 <- tab %>%
# group_by(index, Target) %>%
# slice_max(Percent, n = 4, with_ties = FALSE) %>%
# ungroup()
tab2 <- tab %>%
group_by(index, Target) %>%
filter(Percent > 0.05) %>%
ungroup()
tab_totals <- tab %>%
group_by(index, Target) %>%
summarise(
total = sum(Percent)
)
tab2 <- merge(tab_totals,tab2,
by=c("index","Target"))
left_overs <- tab2 %>% group_by(index,Target) %>%
summarise(
Percent = total - sum(Percent),
) %>% mutate(
index_result = "Other"
) %>% distinct
tab2 <- rbind(tab2 %>% select(-total),left_overs) %>% mutate(
Percent = round(Percent,3)
) %>% filter(
Percent > 0.005
)
tab2$index_result[tab2$index_result=="Gram neg non typeable"] <- "Gram Negative: \n Non-typeable"
tab2$index_result[tab2$index_result=="Gram pos non typeable"] <- "Gram Positive \n Non-typeable"
tab2$index_result[tab2$index_result=="Aerobic non typeable"] <- "Aerobic Growth: \n Non-typeable"
tab2$index_result[tab2$index_result=="E coli"] <- "E. coli"
tab2$index_result[tab2$index_result=="Enterobacter"] <- "Enterobacter spp."
tab2$index_result[tab2$index_result=="Gram neg"] <- "Gram Negative"
tab2$index_result[tab2$index_result=="Gram pos"] <- "Gram Positive"
tab2$index_result[tab2$index_result=="Klebsiella"] <- "Klebsiella spp."
tab2$index_result[tab2$index_result=="Lactococcus"] <- "Lactococcus spp."
tab2$index_result[tab2$index_result=="Nas"] <- "Non-aureus Staph"
tab2$index_result[tab2$index_result=="Other gram pos"] <- "Gram Positive: \n Other"
tab2$index_result[tab2$index_result=="Sslo non typeable"] <- "Strep. or Strep-like \n organism: Other"
tab2$index_result[tab2$index_result=="Contaminated"] <- "Contamination"
#tab2$index_result %>% unique %>% write.csv('index_results.csv')
index_results_order <- read_csv("index_results_order.csv")
tab2 <- merge(tab2,index_results_order,
by="index_result")
tab2 <- tab2 %>% arrange(order)
order <- tab2$index_result %>% unique
tab2$index_result <- fct_relevel(tab2$index_result,"Gram Positive", "Gram Positive \n Non-typeable", "Gram Positive: \n Other",
"Lactococcus spp.", "Non-aureus Staph", "Staph aureus", "Staph chromogenes",
"Staph non typeable", "Strep ag", "Strep dys", "Strep spp", "Strep uberis",
"Strep. or Strep-like \n organism: Other", "Other", "Aerobic Growth: \n Non-typeable",
"Contamination", "Coliform", "Coliform other", "E. coli", "Enterobacter spp.",
"Gram Negative", "Gram Negative: \n Non-typeable", "No growth",
"Klebsiella spp.")
tab2$index_result <- fct_relevel(tab2$index_result,rev)accumast_plot <- ggplot(tab2 %>% filter(index=="Accumast"), aes(x="", y=Percent, fill=index_result)) +
facet_wrap(~ Target,ncol=7) +
geom_bar(width = 2, stat = "identity") +
geom_text(aes(label=paste0(index_result,
" (",
round(Percent*100,1),
"%)")), size=2, position=position_stack(vjust=.5)) +
theme(legend.position = "none",
panel.background = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman"),
axis.ticks.x = element_blank()) +
labs(title = "Accumast") +
scale_y_continuous(labels = scales::percent, breaks = seq(0, 1, 0.2), limits = c(0, NA)) +
xlab("") + ylab("Percentage of cases")
accumast_plot
ggsave("accumast.tiff", width = 8, height = 4, dpi = 200)ggplot(tab2 %>% filter(index=="Biplate"), aes(x="", y=Percent, fill=index_result)) +
facet_wrap(~ Target,ncol=7) +
geom_bar(width = 2, stat = "identity") +
geom_text(aes(label=paste0(index_result,
" (",
round(Percent*100,1),
"%)")), size=2, position=position_stack(vjust=.5)) +
theme(legend.position = "none",
panel.background = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman"),
axis.ticks.x = element_blank()) +
labs(title = "Biplate") +
scale_y_continuous(labels = scales::percent, breaks = seq(0, 1, 0.2), limits = c(0, NA)) +
xlab("") + ylab("Percentage of cases")
ggsave("biplate.tiff", width = 8, height = 4, dpi = 200)ggplot(tab2 %>% filter(index=="Check Up"), aes(x="", y=Percent, fill=index_result)) +
facet_wrap(~ Target,ncol=7) +
geom_bar(width = 2, stat = "identity") +
geom_text(aes(label=paste0(index_result,
" (",
round(Percent*100,1),
"%)")), size=2, position=position_stack(vjust=.5)) +
theme(legend.position = "none",
panel.background = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman"),
axis.ticks.x = element_blank()) +
labs(title = "Check Up") +
scale_y_continuous(labels = scales::percent, breaks = seq(0, 1, 0.2), limits = c(0, NA)) +
xlab("") + ylab("Percentage of cases")
ggsave("checkup.tiff", width = 8, height = 4, dpi = 200)ggplot(tab2 %>% filter(index=="Mastatest"), aes(x="", y=Percent, fill=index_result)) +
facet_wrap(~ Target,ncol=7) +
geom_bar(width = 2, stat = "identity") +
geom_text(aes(label=paste0(index_result,
" (",
round(Percent*100,1),
"%)")), size=2, position=position_stack(vjust=.5)) +
theme(legend.position = "none",
panel.background = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman"),
axis.ticks.x = element_blank()) +
labs(title = "Mastatest") +
scale_y_continuous(labels = scales::percent, breaks = seq(0, 1, 0.2), limits = c(0, NA)) +
xlab("") + ylab("Percentage of cases")
ggsave("Mastatest.tiff", width = 8, height = 4, dpi = 200)ggplot(tab2 %>% filter(index=="Petrifilm"), aes(x="", y=Percent, fill=index_result)) +
facet_wrap(~ Target,ncol=7) +
geom_bar(width = 2, stat = "identity") +
geom_text(aes(label=paste0(index_result,
" (",
round(Percent*100,1),
"%)")), size=2, position=position_stack(vjust=.5)) +
theme(legend.position = "none",
panel.background = element_blank(),
panel.border = element_rect(color = "grey", fill = NA),
text = element_text(family = "Times New Roman"),
axis.ticks.x = element_blank()) +
labs(title = "Petrifilm") +
scale_y_continuous(labels = scales::percent) + xlab("") + ylab("Percentage of cases")
ggsave("Petrifilm.tiff", width = 8, height = 4, dpi = 200)pmap(combinations, top5)[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
[[19]]
[[20]]
[[21]]
[[22]]
[[23]]
[[24]]
[[25]]
[[26]]
[[27]]
[[28]]
[[29]]
[[30]]
[[31]]
[[32]]
[[33]]
[[34]]
[[35]]